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Abstract 

This  study  improved  the  numerical  optimization  technique  used  to  compute  mixed 
H2  / Hoo  controllers.  The  run  time  to  generate  the  full  curve  of  an  H2 /Hoo  design  was  reduced 
by  a  factor  of  almost  200  for  a  specific  example.  This  was  accomplished  by  improving  gradient 
calculations,  upgrading  the  Hoo  norm  calculation  to  capture  multiple  peaks  in  the  singular  value 
plot  and  using  routines  that  minimized  processing  time.  The  optimization  routine  used  was  the 
MATLAB  ’CONSTR’  Sequential  Quadratic  Programming  (SQP)  routine.  This  SQP  routine 
was  compared  to  a  FORTRAN  SQP  numerical  optimizer,  IMSL/DNCONG.  The  method 
was  applied  to  an  F-16  SISO  example  and  a  missile  MIMO  example.  The  mixed  H2IH00 
controller  for  both  these  examples  are  compared  to  previous  controller  designs  created  using 
the  LQG/LTR  technique. 
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IMPROVED  NUMERICAL  SOLUTION  OF  MIXED  H2/H00 
OPTIMIZATION  WITH  APPLICATIONS 


/.  Introduction 

1.1  Objectives 

This  study  will  investigate  improving  computation  of  the  mixed  H2IH0Q  controller 
using  different  optimization  techniques.  The  following  objectives  will  be  addressed: 

1.  Reduce  excessive  run  time  and  improve  reliability  of  the  optimization  code. 

2.  Derive  solutions  for  an  F-16  SISO  example. 

3.  Derive  solutions  for  a  missile  MIMO  example. 

4.  Compare  MATLAB  optimization  routine  to  a  FORTRAN  IMSL/DNCONG  routine. 

1.2  Justification 

Controller  designs  for  aerospace  systems,  which  are  highly  nonlinear,  are  often  based 
on  linear,  time-invariant  models.  The  primary  purpose  of  a  controller  is  to  provide  stability 
and  performance  in  the  presence  of  disturbances  and  uncertainties  (i.e.,  perturbations  or 
unmodeled  dynamics).  The  challenge  to  the  control  designer  is  to  develop  a  linear  solution 
for  these  nonlinear  systems.  A  multi-objective  controller  can  be  used  to  provide  both  stability 
and  performance.  One  example  of  a  multi-objective  controller  is  the  H2IH00  controller.  The 
H2  part  of  the  controller  provides  stability  and  performance  in  the  face  of  noises  acting  on 
the  system.  The  Hoo  part  of  the  design  is  used  to  make  the  system  behave  as  closely  as 
possible  to  some  ideal  linear  model  (model  matching)  and  provides  for  unmodeled  dynamics 
and  uncertainties  in  the  model.  The  goal  then  is  to  find  a  controller  which  minimizes  the  effect 
of  disturbances  while  providing  an  acceptable  level  of  model  matching  and  robustness. 
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The  process  of  computing  an  optimal  H2IH00  controller  is  iterative.  It  is  done  by 
minimizing  the  two-norm  of  one  system  transfer  function  while  constraining  another’s  infinity- 
norm  to  be  less  than  or  equal  to  some  given  value.  Many  different  values  of  the  infinity-norm 
must  be  selected.  The  designer  can  then  make  tradeoffs  in  controller  performance  between 
noise  rejection,  tracking  and  uncertainties  in  the  model.  The  desire  is  to  provide  a  range 
of  controllers  between  the  optimal  H2  and  optimal  ifoo  controllers.  With  the  capability  to 
develop  this  range  of  controllers,  the  designer  can  decide  what  level  of  noise  rejection  versus 
what  level  of  robustness  to  uncertainty  is  required  for  the  system. 

Over  the  past  five  years  several  optimization  routines  have  been  used  to  improve  the 
process  of  generating  a  mixed  controller.  The  first  numerical  optimization  approach  to  find 
a  mixed  H2IH00  controller  was  developed  by  Ridgely  using  a  Davidon-Fletcher-Powell 
(DFP)  routine  [Rid91].  This  routine  used  variable  metrics,  where  information  about  previous 
iterations  are  stored  and  used  to  define  a  better  search  direction  for  the  next  iteration.  This 
took  several  weeks  to  generate  a  set  of  mixed  controllers  when  running  on  a  SPARC  10  SUN 
workstation.  There  were  also  examples  where  the  numerics  caused  such  a  problem  that  the 
mixed  H2IHC0  controller  could  not  be  calculated  for  various  infinity-norm  constraints.  Since 
then,  other  optimization  routines  such  as  penalty  methods  and  conjugate  gradient  routines 
have  been  used  on  this  controller  design  development.  Currently,  the  MATLAB  version  of 
the  Sequential  Quadratic  Programming  (SQP)  method  is  used.  This  method  obtains  a  search 
direction  by  solving  a  quadratic  objective  function  with  linear  constraints.  The  SQP  method 
has  cut  the  run  time  by  a  factor  of  10  from  the  DFP  routine. 

There  were  cases  where  the  infinity-norm  gradients  could  not  be  accurately  calculated 
for  use  in  the  SQP  optimization  routine.  One  of  the  causes  of  this  problem  was  found  to  be 
multiple  peaks  in  the  singular  value  plot  of  the  system.  If  the  infinity-norm  gradients  can 
not  be  accurately  calculated,  the  full  range  of  controllers  can  not  be  designed.  Without  this 
full  range  of  designs,  the  engineer  is  limited  in  the  choice  of  controllers.  Algorithms  have 
been  developed  in  this  work  that  take  into  account  the  multiple  peak  problem  and  allow  the 
calculation  of  the  mixed  controller. 
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There  are  two  specific  control  design  examples  for  which  the  mixed  controllers  could 
not  be  calculated  at  a  full  range  of  infinity-norm  values.  These  include  an  F-16  Single-Input- 
Single-Output  (SISO)  normal  acceleration  command  following  controller  and  a  controller 
for  a  tail/fin  guided  missile  that  is  a  Multiple-Input-Multiple-Output  (MIMO)  set-up.  The 
algorithm  developed  to  improve  the  accuracy  of  the  Hoo  gradient  will  be  used  on  these  two 
examples. 

The  new  algorithm  is  much  faster  and  more  reliable  than  the  previous  version.  This 
provides  the  engineer  with  a  more  practical  and  complete  method  for  designing  controllers  for 
uncertain  (possibly  nonlinear)  systems  operating  in  the  face  of  disturbances. 

1 3  Literature  Review 

The  initial  approach  to  the  mixed  objective  problem  was  developed  by  Bernstein  and 
Haddad  [BH89].  Their  approach  was  limited  in  that  it  minimized  an  upper  bound  to  the 
two-norm  rather  than  the  two-norm  itself.  It  also  lacked  the  ability  to  define  the  order  of  the 
controller  and  to  freely  define  the  two-norm  and  infinity-norm  transfer  functions.  Work  done 
by  Zhou  [ZDB89],  [ZDGB90]  presents  an  approach  where  different  exogenous  inputs  were 
allowed,  but  only  one  controller  output.  Yeh  [YBC90]  showed  this  to  be  the  dual  of  [BH89]. 

The  mixed  controller  design  problem  of  [BH89]  was  transformed  into  a  convex  nu¬ 
merical  optimization  program  by  Khargonekar  and  Rotea  [KR91].  This  approach  is  an 
effective  alternative  to  solving  coupled  nonlinear  equations  required  by  the  method  developed 
by  [BH89].  Several  software  routines  are  currently  available  to  solve  convex  optimization 
problems.  Rotea  and  Khargonekar  were  the  first  to  develop  a  true  mixed  f/^z/^oo  solution, 
valid  only  for  a  special  case.  They  did  not  use  an  upper  bound  to  the  two-norm  but  rather  the 
two-norm  itself.  They  also  allowed  for  two  sets  of  exogenous  inputs  and  controlled  outputs. 
The  one  restriction  of  their  development  was  full  state  availability.  Ridgely  developed  a 
numerical  solution  to  nine  necessary  conditions  for  an  optimal  mixed  H2I Hoo  controller.  He 
solved  an  unconstrained  optimization  problem  for  the  controller  state  space  matrices  using 
the  DFP  algorithm  [Rid91].  His  work  also  lifted  the  restriction  of  full  state  availability  and 
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developed  the  conditions  that  an  output  feedback  H2IH00  controller  must  satisfy.  The  most 
recent  contribution  to  the  solution  of  this  problem  was  by  Walker  [Wal94],  who  significantly 
reduced  the  time  required  to  generate  the  H2IH00  controller.  He  reformed  the  mixed  H2IH00 
optimization  design  process  into  a  constrained  optimization  problem  with  an  Hoo  constraint. 
He  then  applied  the  MATLAB  SQP  method  of  optimization  to  the  constrained  nonlinear 
problem. 

Walker  identified  the  problem  of  the  multiple  peaks  in  the  singular  value  plot  as  the 
reason  some  problems  do  not  converge  to  a  mixed  controller  [Wal94].  The  infinity-norm  is 
defined  by  doing  a  frequency  search  of  the  magnitude  of  the  singular  value  plot.  The  point 
where  the  maximum  occurs  is  the  infinity-norm.  The  difficulty  in  calculating  the  infinity- 
norm  gradients  is  created  when  there  is  more  than  one  frequency  at  which  the  maximum 
occurs.  The  optimization  routine  does  not  know  at  which  frequency  the  norm  should  be 
differentiated.  It  will  continually  shift  back  and  forth  from  one  frequency  to  another,  and  thus 
prevent  convergence  to  a  mixed  solution. 

Several  types  of  examples  have  been  used  to  demonstrate  the  mixed  controller  design 
approach.  Luke  [Luk93]  used  an  example  of  an  F-16  normal  acceleration  command  following 
model.  His  results  showed  that  controllers  designed  using  mixed  H2IH00  provided  excellent 
performance  and  robustness  at  orders  less  than  the  order  of  the  conventional  LQGA-TR  control 
design.  He  also  addressed  the  problem  of  a  MIMO  tail-fin  controlled  missile.  He  used  the 
MATLAB  SQP  routine,  but  the  problem  became  numerically  ill  conditioned  with  no  useful 
mixed  controllers  calculated. 

1.4  Approach 

Reduction  of  the  run  time  for  development  of  the  mixed  H2/H00  controller  was  at¬ 
tempted  first  by  trying  other  optimization  routines.  The  original  software  developed  by 
Walker  was  run  using  the  MATLAB  optimization  toolbox  [Gra93].  MATLAB  has  one  built-in 
optimization  routine,  which  is  ’CONSTR’.  There  are  several  other  routines  written  in  FOR¬ 
TRAN  that  are  numerically  well  conditioned  for  this  type  of  problem  and  that  can  provide 
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faster  run  times.  MATLAB  can  be  linked  to  FORTRAN  to  provide  the  capability  of  running 
FORTRAN  optimization  routines.  One  of  the  FORTRAN  routines  that  will  be  tested  is  the 
DNCONG  routine  [IMS89].  DNCONG  is  a  double  precision  routine  that  uses  the  SQP  al¬ 
gorithm  with  user  supplied  gradients.  This  will  allow  a  comparison  between  the  line  search 
and  Hessian  update  methods  used  in  MATLAB  SQP  with  the  the  methods  used  in  FORTRAN 
IMSL  SQP  routine. 

Other  things  were  also  done  to  optimize  the  run  time.  First,  most  of  the  time  used  to 
calculate  the  mixed  controller  is  spent  doing  the  frequency  search  to  find  the  infinity-norm. 
This  search  can  be  optimized  by  doing  one  general  search  to  find  the  peaks  of  the  problem, 
and  then  reducing  the  search  down  to  a  smaller  range  of  frequencies  about  those  peaks.  This 
should  not  only  improve  run  times  but  should  also  improve  the  accuracy  of  the  infinity-norm 
and  its  gradients.  Second,  an  active  set  strategy  can  be  applied  to  only  calculate  the  derivatives 
of  the  constraint  when  the  program  determines  that  the  constraint  is  active.  This  is  estimated 
to  be  where  the  second  largest  amount  of  computer  time  is  spent.  Third,  the  run  time  may 
also  be  improved  by  defining  less  storage  space  and  having  fewer  function  calls.  Finally, 
normalizing  the  objective  and  constraints  should  also  improve  the  efficiency  of  the  code. 

The  most  important  improvement  to  the  optimization  code  was  the  correction  of  the 
gradients  for  the  objectives  and  constraints.  Modifications  were  made  to  the  gradients  to 
take  into  account  the  subdiagonal  terms  of  the  controller  state  space  matrices  once  it  was 
put  into  modal  form.  Corrections  also  included  a  change  to  the  two-norm  gradient  when  the 
mixed  controller  was  not  at  the  H2  optimal  design.  Modifications  were  also  made  to  the 
stability  constraint  gradient  so  that  it  could  be  calculated  analytically  instead  of  using  finite 
differencing. 

The  problem  of  multiple  peaks  when  calculating  the  Hoa  constraint  was  also  addressed. 
The  solution  was  to  incorporate  a  routine  that  defines  where  the  peaks  occur.  This  can  be  done 
by  checking  the  slopes  of  the  singular  value  plot.  Peaks  can  be  identified  wherever  the  slope 
of  the  singular  value  plot  is  zero,  and  the  slope  changes  fi-om  positive  to  negative  about  that 
point.  An  active  set  strategy  can  then  be  applied  to  constrain  the  problem  at  each  point  where 
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a  peak  is  identified.  Using  separate  constraints  for  each  peak  will  prevent  the  optimizer  from 
jumping  back  and  forth  between  peaks  and  causing  the  problem  to  not  converge. 

Other  improvements  were  made  to  the  code  to  reduce  the  run  time,  prevent  non¬ 
stabilizing  controller  designs,  and  provide  the  capability  to  step  along  the  mixed  H2IH00 
curve  automatically.  A  function  in  MATLAB  called  FMIN  was  used  to  find  maximum  points 
in  the  singular  value  plot  once  a  peak  had  been  defined.  This  was  much  more  efficient  than  the 
several  hundred  point  search  over  the  entire  frequency  range.  Also,  to  prevent  non-stabilizing 
solutions  a  penalty  function  was  added  to  the  objective  function  when  the  stability  constraint 
was  violated.  This  proved  to  be  an  effective  means  to  prevent  non-stabilizing  designs.  An 
autostep  function  was  also  developed  that  allows  the  designer  to  set  a  desired  percentage 
change  in  two-norm  and  infinity-norm  and  the  function  automatically  steps  along  the  curve. 
This  reduced  the  run  time  and  the  required  interaction  by  the  designer  with  the  code. 

Finally,  the  inability  to  produce  a  mixed  solution  for  some  specific  examples  will  be 
addressed.  The  algorithms  developed  for  the  multiple  peak  problem  will  be  used  along  with 
other  improvements  in  the  code  to  provide  a  full  range  of  solutions  to  these  cases.  The  SISO 
F-16  acceleration  command  following  problem  and  the  MIMO  tail/fin  controlled  missile 
problem  have  been  identified  as  cases  with  multiple  peaks.  Application  of  the  algorithm 
should  provide  improvements  in  the  capability  of  calculating  the  mixed  controllers  and  reduce 
the  time  required  to  find  them. 
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11.  Background 


This  chapter  discusses  the  theory  used  to  develop  the  mixed  H2IH00  controller  along 
with  background  of  numerical  optimization.  The  discussion  on  the  mixed  H2IH00  controller 
is  broken  into  three  parts.  The  first  part  begins  with  a  discussion  on  how  the  H2  controller  is 
developed  and  then  the  background  of  the  Hoo  controller  follows.  After  the  separate  design 
methods  for  H2  and  if 00  have  been  presented,  an  explanation  of  how  these  two  methods  are 
combined  to  form  a  mixed  design  is  given.  The  numerical  optimization  section  gives  a  brief 
overview  of  optimization,  presents  the  history  of  the  numerical  optimization  of  the  mixed 
if2/ifoo  controller  and  provides  a  discussion  of  the  SQP  optimization  routine.  The  numerical 
optimization  section  also  contains  the  methods  used  to  calculate  derivatives  of  the  H2  and 
Hoo  norms. 

2.1  Mixed  H2 / if 00  Controller  Design 

2.1.1  H2  Optimization.  The  H2  method  is  a  well  defined  process  to  design 
internally  stabilizing  controllers  for  linear  dynamic  systems.  The  objective  of  an  H2  design  is 
to  synthesize  a  controller  which  minimizes  the  energy  of  a  system  output  to  a  zero-mean  white 
Gaussian  noise  input.  The  H2  method  is  a  generalization  of  the  standard  Linear  Quadratic 
Gaussian  (LQG)  problem. 

The  H2  optimization  process  can  be  explained  using  the  diagram  shown  in  Figure  2.1. 
The  linear  time-invariant  model  of  the  system  dynamics  is  defined  by  P.  The  compensator 
developed  through  H2  design  is  represented  by  K.  Inputs  to  the  system  are  zero-mean.  Unit 
intensity  White  Gaussian  Noise  (UWGN)  shown  as  w.  The  measured  values  coming  out  of 
the  system  are  shown  as  y  and  are  input  into  the  compensator.  The  control  signal  generated 
by  the  compensator  is  labeled  u  and  is  fed  back  into  the  system.  The  controlled  output  of  the 
system  is  represented  as  z.  The  goal  of  this  method  is  to  design  a  K  that  minimizes  the  energy 
of  the  output  2:,  given  a  (UWGN)  w. 
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The  plant  represented  as  P  in  the  diagram  can  be  partitioned  into  its  separate  input  and 
output  transfer  functions 


P  = 


P  P 

^  zw  zu 


p  p 

yw  yu 


(2.1) 


where 

Z  —  P^^yjIV  “f*  Pzu'^ 

y  —  Pyw  "H  ^yu  ^ 

A  state  space  realization  of  P  can  be  written  as 


(2.2) 


X2  =  A2X2  +  B^w  +  Bu2U 

Z  =  CzX2  +  Dzw'W  +  DzuU  (2.3) 

y  —  ^1/2^2  Byw"^  “I”  Byu'^ 


Minimizing  the  two-norm  of  the  closed-loop  system  is  the  same  as  minimizing  the  energy  of 
the  output  of  the  system,  given  a  (UWGN)  input. 


a 


inf,,  |1;^1|2  foriu  6  UWGN 

K(s)  Stabilizing 
K(s)}ubilizing  ^ 

1 ■  ll-P^w  "I"  PzuK{I  —  PyuK)  -Pj/wlh 

K  (s)biaoiliztng 


The  following  assumptions  are  standard  when  designing  an  H2  controller: 


(2.4) 

(2.5) 

(2.6) 
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(i)  =  0 

(ii)  Dyu  =  0 

(iii)  ( ^2,  Bu2  )  is  stabilizable  and  {Cy^ ,  A2)  is  detectable 

(iv)  rank  and  rank 


A2  -  joji 

Bu2 

(V) 

has  full  column  rank  for  all  uj 

c. 

D,u 

(Vi) 

A2  -  juji 

Bn, 

has  full  row  rank  for  all  uj 

Cy2 

Dynj 

Assumption  (i)  is  required  so  that  the  closed  loop  two-norm  of  the  system  is  finite. 
Condition  (ii)  makes  the  development  easier  but  can  be  removed  completely.  Requirement 

(iii)  is  required  for  existence  of  a  stabilizing  solution.  This  guarantees  that  the  set  of  stabilizing 
controllers  is  not  empty.  Assumption  (iv)  is  needed  to  prevent  the  use  of  inifinite  control  power 
and  avoids  the  singular  control  problem.  Finally,  requirements  (v)  and  (vi)  are  required  to 
ensure  the  existence  of  solutions  to  the  two  Algebraic  Ricatti  Equations  (ARE)  required  in  the 
H2  solution.  These  AREs  are  defined  in  (2.14)  and  (2.16). 

All  H2  controllers  can  be  represented  by  a  Lower  Fractional  Transformation  (LFT) 
of  a  transfer  function  J  and  a  freedom  parameter  Q  as  shown  in  Figure  2.2.  The  freedom 
parameter  is  used  to  identify  suboptimal  H2  controllers.  These  suboptimal  controllers  may  be 
useful  when  used  for  computing  designs  that  trade  off  H2  performance  for  if 00  performance. 
The  mixed  H2/H00  optimization  process  will  usually  start  at  the  H2  optimal  solution,  where 
Q=0.  Rescaling  y  and  u  so  that  condition  (iv)  is  strengthened  to 

(iv) /  =  /  and  Dy^D^^  =  I 

results  in  the  following  H2  optimal  controller  parameterization: 
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Figure  2.2  system  with  parameterized  controller 


where 


Juy 

J iir 

Aj 

Kf 

Kfi 

As)  = 

— 

-Kc 

0 

I 

Jvy 

Jyr 

Kci 

I 

0 

(2.7) 

(2.8) 


Aj  =  A  —  KjCy^  —  Bu2Kc 

(2.9) 

a;  =  -SJ,A-2  +  PJ.C. 

(2.10) 

K,  =  Y2CI  +  B.DI, 

(2.11) 

II 

1 

(0 

(2.12) 

Kji  =  Bu2 

(2.13) 

X2  and  Y2  are  the  stabilizing  solutions  to  the  AREs 

{A  -  Bu,DlC,fX2  +  X2{A  -  -  X2Bu,BIX2  +  CjC,  =  0  (2.14) 


where 

C,  =  {I-  (2.15) 
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and 


(A  -  B,DI,C„)Y2  +  Y,{A  -  -  Y,CIC„Y2  +  B.Bl  =  0  (2.16) 

where 

=  B^I  -  Dl^Dy^,)  (2.17) 

If  <5  =  0,  the  controller  is  unique  and  optimal.  If  Q  is  nonzero,  the  set  of  all  controllers  such 
that  ||r2,„||  <  a  is  parameterized  by  any  Q  6  RH^  such  that 

\\Q\\l<c^-a^  (2.18) 

Another  approach  to  minimize  the  two-norm  of  the  system  is  through  the  use  of  the 
definition  of  the  controllability  grammian  and  a  Lyapunov  equation.  Assume  the  state  space 
closed-loop  A  matrix  for  the  H2  problem  is  stable.  Define  the  controllability  grammian  as 

Lc  =  /  (2.19) 

Jo 

where  Lc  is  the  unique,  real,  symmetric  positive  semi-definite  solution  to  the  Lyapunov 
equation 

0  =  ALc  +  LcA'^  +  BB^  (2.20) 


The  two-norm  is  then  defined  as 


\\G{M\\l  =  tr[CLcC^]  (2.21) 

This  method  was  used  for  the  mixed  H2IH00  optimization  code.  The  Lyapunov 
equation  is  solved  using  the  MATLAB  command  ’GRAM’.  The  objective  function  used  for 
the  optimization  process  is  then  the  square  of  the  two-norm,  with  the  Lyapunov  equation 
appended  as  a  constraint. 
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Figure  2.3  Perturbed  stability  diagram 


2.1.2  Hoo  Optimization.  The  objective  of  an  Hoo  controller  is  to  minimize  the 
energy  of  the  output  to  the  worst  possible  bounded  energy  input.  This  is  the  same  as  minimizing 
the  infinity-norm  of  the  transfer  function.  This  type  of  design  is  robust  to  uncertainties  in 
the  plant.  Let  the  uncertainty  be  called  A.  Figure  2.3  shows  how  the  uncertainty  would  be 
introduced  to  the  system.  The  question  of  how  large  a  A  can  be  introduced  and  the  system 
remain  stable  is  answered  by  use  of  the  small  gain  theorem  [Zam89].  For  the  closed  loop 
system  to  be  stable,  Equation  (2.22)  must  be  satisfied. 

Theorem  1  (Small  Gain  Theorem)  Let  Ted  €  Hoo.  Assume  A  e  Hoo  is  connected  from  e 
to  d  as  shown  in  Figure  23.  Then  the  closed-loop  system  is  internally  stable  if 


||red(5)A(5)||oo  <  ||re<i(5)|U||A(5)||oo  <  1  (2.22) 


Proof:  [Zam89] 

Thus,  by  minimizing  ||Ted||oo,  the  allowable  ||A||oo  that  cannot  destabilize  the  system  is 
maximized.  Hoo  optimization  is  given  by 


inf  sup  ||e||2  =  inf  UTedHoo  =  7  (2.23) 

^adm  ||ci||2<l  ^adm 
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Figure  2.4  H^o  Feedback  System 


where  the  infinity-norm  of  Ted  is 


\\Ted\\oo  =  SUpcrfTed] 

<JJ 


(2.24) 


and  a  denotes  the  maximum  singular  value.  The  minimum  achievable  infinity-norm,  as 
indicated  by  Equation  (2.23),  is  designated  7. 

The  design  diagram  for  Hoc  is  the  same  as  the  figure  for  H2  except  for  the  definition  of 
the  inputs  and  outputs.  Where  the  input  to  H2  was  a  white  noise  labeled  w,  the  input  for  H^o 
is  a  bounded  energy  exogenous  input  labeled  d.  The  output  for  H2  was  a  controlled  output  z, 
where  the  controlled  output  for  Hoo  is  labeled  e.  This  diagram  is  shown  in  Figure  2.4. 


by 


or 


The  plant  P  in  Figure  2.4  can  be  partitioned  into  its  separate  transfer  functions  as  given 


P  = 


Ped  Peu 


(2.25) 


e  —  Pedd  -j-  P eu'^ 
y  —  Pydd  -|-  Pyu^ 


(2.26) 
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A  state  space  realization  of  the  plant  P  is  defined  by 


^oo 

8 

H 

8 

Ii 

+ 

Bdd 

+ 

By^^U 

e 

—  ^e^oo 

+ 

+ 

Day'll 

y 

+ 

Bydd 

+ 

Byu  U 

The  following  assumptions  are  made  for  the  Hoo  design: 


(2.27) 


(i)  Ded  =  0 


(ii)  Dyu  =  0 


(iii)  (Aoo,  Bu^)  is  stabilizable  and  iCy^,Aoo)  is  detectable 


(iv)  Dj^Deu  =  /  and  DydDy^  =  I  (strengthened  through  scaling) 


(V) 


-^OO  Buca 

Cg  Dgy, 


has  full  column  rank  for  all  uj 


(Vi) 


Aoc  -  joji  Bd 

^y<x  Byd 


has  full  row  rank  for  all  u) 


Conditions  (i)  and  (ii)  are  not  required  for  a  solution  to  exist,  but  reduce  the  complexity  of 
the  solution.  Condition  (iii)  is  necessary  for  the  existence  of  a  stabilizing  controller.  Condition 
(iv)  is  a  requirement  that  all  control  usage  is  penalized  and  no  perfect  measurements  are 
allowed.  Finally,  conditions  (v)  and  (vi)  guarantee  there  is  a  solution  to  the  Riccati  equations 
for  the  Hoo  norm  calculation. 


The  Hoo  controller  is  found  by  an  iterative  process.  The  Hoo  optimal  controller  found 
by  this  process  is  not  unique.  The  family  of  all  Hoo  compensators  such  that 


ll^ecilloo  <  7 


(2.28) 


is  given  by  (see  Figure  2.5) 


i^(s)  =  F,[J(^),Q(.)] 


(2.29) 
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where 


Figure  2.5  i^oo  system  with  parameterized  controller 


J uy 

Jut 

Aj 

Kf 

Kfi 

J{s)  = 

— 

-Ko 

0 

I 

J vy 

Jvr 

Kci 

I 

0 

K,  =  {Bu^X^  +  DlC,){I-r"Y^Xoo)-^ 

Kf  =  Y^Cj^+B,Dl, 

Kci  =  -{r^Dy,BjX^  +  CyJ{I-r"Y^X^)-^ 

Kf,  =  r"Y^CjD,u+B^^ 

The  matrices  Xoo  and  are  the  solutions  to  the  AREs 


{A  -  BuDlC^fX^  +  Xoo{A  -  Bu^DlC,) 

+Xooir^B,Bj  -  B^^BljX^  +  CjCe  =  0 


(2.30) 


(2.31) 

(2.32) 

(2.33) 

(2.34) 

(2.35) 


(2.36) 


\ 
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where 

C,  =  {I-DeuDl)Ce  (2.37) 

and 

(A  -  +  Y^{A  - 

+I"„(7“"CjC,  -  Cl^C,^)Y^  +  BiBj  =  0  (2.38) 


where 


Bd  =  Bd{I  -  Dy^Dyd) 


(2.39) 


Finally,  Q  can  be  chosen  to  be  any  Q  G  Boo  such  that 

WQWoo  <  7  (2.40) 

The  above  parameterization  of  a  controller  K  is  valid  if  and  only  if  the  following  three 
conditions  hold; 

(i)  Hx  €  dom{Ric)  with  Xoo  =  Ric{Hx)  >  0 

(ii)  By  €  dom{Ric)  with  Y^o  =  Ric{BY)  >  0 

(iii)  /9(FooXoo)  <  7^ 

The  process  of  finding  a  controller  which  results  in  a  closed-loop  infinity-norm  close 
to  7  is  iterative.  To  find  a  controller,  select  7  and  check  to  see  if  it  meets  all  three  conditions 
above.  If  any  condition  fails,  increase  7  and  repeat  the  check.  If  all  three  conditions  are  met, 
decrease  7.  The  desired  7  can  be  found  by  this  process  to  any  desired  accuracy. 

2.13  Mixed  B2  /if 00  Optimization.  A  method  was  derived  to  achieve  the  robusmess 
of  if 00  control  with  the  noise  rejection  properties  of  B2  by  Ridgely  [Rid91].  The  mixed 
B2  / if 00  diagram  is  shown  in  Figure  2.6.  Using  this  method,  weights  can  be  added  to  both  the 
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d 

w 

u 


K 

Figure  2.6  General  mixed  H2 / //^oo  optimization  problem 

H2  and  Hoo  set-up.  The  final  goal  of  this  mixed  design  process  is  to  develop  an  admissible 
controller  that  satisfies  the  following: 


inf  ,  ||r2^||2  subject  to  ||Te<i||oo  <  7 

h  staothztng 


A  general  form  of  the  mixed  system  can  be  written  as 


e 

Bed  Pew  Peu 

d 

2 

= 

Pzd  P zw  P zu 

W 

y 

Pyd  Pyw  Pyu 

u 

The  state  space  form  of  the  system  can  be  written  as 

X  =  Ax  -f  Bad  -h  ByjW  BuU 


e  —  C qX  -(”  B^^d 
z  —  C zX  -j-  B ^^d  -{“  B -{-  B 

y  -  CyX  -j-  Bydd  “1“  By-ijuXO  “j”  B ynXl 


\ 


(2.41) 


(2.42) 

(2.43) 

(2.44) 

(2.45) 
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The  individual  H2  and  H^o  problems  can  be  approached  separately,  as  described  in 
Subsections  2.1.1  and  2.1.2.  The  goal  for  H2  is  to  find  an  admissible  controller  that  minimizes 
the  two-norm  of  while  the  goal  for  Hoc,  is  to  find  an  admissible  controller  for  Ted  that  has 
an  infinity-norm  less  than  or  equal  to  some  specified  7.  The  closed-loop  transfer  functions 
Ted  and  T^w  for  the  mixed  problem  are  defined  as 


=  C,isI-A2)-^B^  +  V,^  (2.46) 

Ted  =  Ce{sI-Aoo)-^Bd  +  Ved  (2.47) 

The  state  space  formulation  of  the  mixed  i/2/^00  controller  is  defined  as: 


Xe  —  AqXo  -j-  HeU 
U  =  CcXe 

The  individual  closed  loop  state  space  matrices  can  be  represented  by 


X2 


Xno  - 


— 


B„ 


X2 

Xc 

X<x> 

Xr 


A2  Bu2Cc 
BcCy2  Ac 


Aoo  By^ooCc 


BcT'yoo  Ac 


B„ 


BcBy^ 


(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 

(2.53) 

(2.54) 
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Bd  = 

C,  = 

Ce  = 
= 

'I^ed  = 


Bd 

BcDyd 


Cz  DzuCc 


Ce  DeuCc 


0 

Bed 


(2.55) 

(2.56) 

(2.57) 

(2.58) 

(2.59) 


The  assumptions  for  the  mixed  problem  can  simply  be  a  combination  of  the  assumptions 
for  the  separate  H2  and  Hoo  problems.  However,  Walker  [Wal94]  has  reduced  the  required 
assumptions  to  the  following: 


(i)  Dzn.  =  0 

(ii)  Dyu  =  0 

(iii)  (^2,  Bu2)  stabilizable,  {Cy^^A^)  detectable 


(iv)  DJ^Dzu  full  rank,  Dy^Dy^  full  rank 


(V) 

(Vi) 


A2  -  jul  Bu2 
Cz  Dzu 

A2  -  jU}I  Byj 
r  n 

'^3/2  -^yu! 


has  full  column  rank  for  all  u; 


has  full  row  rank  for  all  u 


These  are  simply  the  assumptions  for  the  H2  problem.  Since  (i-vi)  guarantee  a  strictly 
proper  admissible  controller.  Walker  found  that  the  H^o  assumptions  could  be  relaxed.  Ded 
no  longer  must  be  zero  and  no  assumptions  need  to  be  made  on  the  ranks  of  D^u  and  Dyd  (i.e., 
singular  and  non-strictly  proper  Hoo  constraints  are  allowed). 

The  following  definitions  are  made  to  simplify  the  development: 


inf 

K  admissible 


(2.60) 
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lU 


Kr 


a  := 

opt  * 

7  := 


7 

a* 


K  -hi 

hadmtsstble 

the  unique  K{s)  that  makes  ||T^u,||2  =  a 
llTedlloo  when  K{s)  = 

a  solution  to  the  H2IH00  problem  for  some  7  >  7 

ll^erflloo  when  —  Kmix 
||7'zu;||2  when  K(^S^  —  Kmix 


(2.61) 

(2.62) 

(2.63) 

(2.64) 

(2.65) 

(2.66) 


where  "admissible"  in  (2.60)  and  (2.61)  means  stabilizing  and  at  a  fixed  order.  The  conditions 
of  the  mixed  H2IH00  problem  can  now  be  stated  as: 

(i)  -42  and  Aoo  are  stable 

(ii)  ||red||oo  <  7  for  some  given  7  >  7 

(iii)  1 1  Tzv,  1 1 2  is  minimized. 

The  following  theorem  defines  the  problem  further: 

Theorem  2  Let  (Ac,  Be,  Cc)  be  given  and  assume  there  exists  a  Qoo  =  >  0  satining 

AeoQoo  +  QooAl  +  {Q^Cj  +  BdVj^)R-\QeoCj  +  +  B^BJ  =  0  (2.67) 

where  R  =  (7^/  —  Ve/Dlj)  >  0.  Then  the  following  are  equivalent: 

(i)  ( Aoo  ,Bd)  is  stabilizable 

(ii)  Aoo  is  stable 

(iii)  A2  is  stable. 


Moreover,  if  the  above  hold  then  the  following  are  true: 


(v)  the  two-norm  of  the  transfer  function  is  given  by 

||r„||5  =  <r[C.QjCjl  =  tTlQ2C%] 
where  Q2  =  ^  ^  solution  to  the  Lyapunov  equation 

•^2*52  +  Q2-^  +  =  0 

(vi)  all  real  symmetric  solutions  Qoo  of  (2.67)  are  positive  semidefinite 

(vii)  there  exists  a  unique  minimal  solution  Qoo  to  (2.67)  in  the  class  of  real  symmetric 
solutions 

(viii)  Qoo  is  the  minimal  solution  of  (2.67)  iff 
Re[\i{Aoo  +  BdV^^R-^Ce  +  QooCjR-K)]  <  0  for  alii 

(ix)  \\Ted\\oo  <  i<)7iff^e]^\i{Aoo  +  BdT>J^R~'^Ce-\-QooCj'R~'^Ce)  <{<)Owhere 
Qoo  is  the  minimal  solution  to  (2.67). 

Proof:  See  [Wal94]  Theorem  (5.2.1) 

The  result  of  this  theorem  is  that,  given  a  controller  which  satisfies  (2.67)  [and  thus 
bounds  IITedlloo],  we  can  determine  the  minimum  a  by  also  requiring  tx[CzQ2Cj\  to  be 
minimized.  The  mixed  problem  can  then  be  rewritten  as:  determine  the  controUer(Ac,  Be  and 
Cc)  that  minimizes 

J{Ao,Bo,Co)  =  ||T,„||2  =  tr[Q2CjC,]  (2.68) 

where  <52  is  the  real  symmetric,  positive  semidefinite  solution  to 

>12^2  +  <32^2  +  =  0  (2.69) 


2-15 


and  such  that  there  exists  a  real,  positive  semi-definite  solution  to 


+  QooAl  +  (QooCj  +  B,Vl)R-'{Q„Cj  +  +  BSl  =  0  (2.70) 

The  minimization  problem  now  looks  like  an  objective  function,  (2.68),  and  two  con¬ 
straint  functions,  (2.69)  and  (2.70).  To  solve  this  minimization  problem  a  Lagrange  multiplier 
approach  is  used.  The  Lagrangian  is 

C  =  tr[Q2C^Cz]  +  tr{[A2Q2  +  Q2-^2  + 

+  tr{[AooQoo  +  QooAl  +  (QooCj  +  Bd'Dl)R-\QooCj  + 

+  B,Bj]y}  (2.71) 

This  approach  was  introduced  by  Ridgely  [Rid91]  and  modified  by  Walker  [Wal94]  to  include 
the  Ded  term  and  allow  for  the  singular  H^o  constraint.  The  first  order  necessary  conditions 
for  the  minimum  of  the  Lagrangian  (where  X  and  3^  are  symmetric  Lagrange  multipliers)  are 


dC 

dAc 

dC 

dB, 


dC 

dC, 


—  ^12^12  +  X2Q2  +  Y^Qab  +  Y2Qb  —  0  (2.72) 

=  +  X2QI2CI  +  X^,V,2  +  X2B,V2  + 

+  Y2QLcJ^  +  Y^^V^b  +  Y2BM  +  {Y^2Qa  +  Y2QI)CJM 

(Y.lQab  +  Y2Qb)CjDlM  =  0  (2.73) 

=  B^^X\Qi2  +  B^^Xi2Q2  +  RuQi'i  +  R2CCQ2  +  B^^YiQab 
Y  B^^Y\2Qb  +  RabQoYlQab  +  R-abQaYuQb  +  RabQabY^Qab 
+  R^bQabY2Qb  +  RbCcQ^bYlQab  +  RbCcQbY^Qab 
+  RbCcQ^bYuQb  +  RbCcQbY2Qb 

+  PliYiQab  +  Yi2Qb)  +  P2{Y^2Qab  +  ^aQi)  =  0 


—  >^2Q2  +  Q2-^2  +  —  0 


dC 
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(2.74) 

(2.75) 


dc 

W2 

dy 

dC 

dQoo 


=  +  ?!:A2  +  cJc^^o 


(2.76) 


—  -^ooQoo  -\- Qoo'^'^  +  (QooCj  +  ^{Qco(^I  ^  ^dVed)^ 

4-  BdBj  =  0  (2.77) 

=  (Aoo+B.V'^^R-^Ce  +  QooCjR-^Cefy 

+  y{Aoo  +  BdVj^R-^C,  +  QooCj  R-^Ce)  =  0  (2.78) 


where 


M  =  R-^DedD'^d 

Pi  =  DlR-^DedBj 


P2  = 

Q2  = 

A’  = 

Qoo  = 

y  = 


dImbJ 
Ql  Qi2 
QL  Q2 

Xt  Xi2 

Xl  X2 

Qa  Qah 
Qjb  Qf> 


Yi  Fi 


12 


Y^l  Y2 


(2.79) 

(2.80) 
(2.81) 

(2.82) 

(2.83) 

(2.84) 

(2.85) 


B^B 


T 

w 


p„ 


BcDyyj 


Bl 


V,  yi2Pj 

B,VT  BcV^BJ 


(2.86) 


+ /)Bj  = 

Bd 

+ 1) 

Bj  D^dBj 

^c^yd 
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(2.87) 


K 

bmbJ 


cl 

CjDl 


Cz  Dzu  Cc 


CjR-^C,  = 


i?l  R\2Cc 

CjRl^  CjR^a 

Cl 

CTDI. 


C  --2'-'c 

-1 


(2.88) 


R 


a  D  e.u  Cc 


Ra  RahC c 


CjRkCc 


(2.89) 


Equations  (2.75)  and  (2.77)  are  the  original  constraint  functions.  One  solution  to  (2.78) 
is  that  y  equals  zero;  the  other  possibility  is  that  (^oo  +  BdV'^^R~'^Ce  +  QooC^ R~^Ce) 
is  neutrally  stable.  The  first  solution  corresponds  to  7  being  off  the  boundary  of  the  Hoc 
constraint.  The  second  solution  corresponds  to  7*  lying  on  the  boundary  of  the  Hcx>  constraint. 
The  off  boundary  condition  of  =  0  results  in  the  Lagrangian  reducing  to 


C  —  tr\Q2C^Cz\  +  ^^{[.42^2  +  ^2^4^  +  (2.90) 

which  is  just  the  H2  Lagrangian;  thus,  for  iV  =  0,  Kmix  =  K2opt-  The  following  theorem 
was  developed  by  Walker  [Wal94]. 

Theorem  3  Assume  Uc  >  122,  where  Uc  is  the  order  of  the  controller  and  n2  is  the  order  of 
the  plant  for  the  H2  problem.  Then  the  following  hold: 

1.  If  ^  <  jj  no  solution  to  the  mixed  H2IH00  problem  exists 

2.  Ifj_  <  7  <  7,  Kmix  is  such  that  7*  =  7 

3.  If  7  >7,  K2^t  is  the  solution  to  the  mixed  H2IH00  problem. 
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Figure  2.7  Typical  H2IH00  curve 

In  summary,  the  solution  to  the  mixed  H2  /  Hoo  problem  will  lie  on  the  Hoo  constraint 
boundary  whenever  the  order  of  the  controller  is  chosen  greater  than  or  equal  to  that  of  the 
H2  problem.  Theorem  3  implies  that,  for  7  <  7  <7,  each  point  on  the  mixed  H2I Hoo 
curve  occurs  on  the  boundary  of  the  Hoo  constraint.  Therefore,  between  7  and  7  the  value 
of  a*  is  a  smooth  monotonically  decreasing  function  as  shown  in  Figure  2.7.  This  type  of 
curve  assumes  the  optimization  process  has  found  the  global  minimum  of  the  mixed  H2I  Hoo 
problem  or  that  the  local  minimum  does  not  change  between  7  steps.  Currently,  the  mixed 
problem  can  not  be  solved  analytically,  and  must  be  solved  through  a  numerical  approach. 
The  numerical  approach  is  addressed  in  Section  2.2. 

2.2  Numerical  Optimization 

2.2.1  Optimization  Background.  This  section  is  a  general  overview  of  optimization 
theory.  For  details  the  reader  is  referred  to  [Van84].  The  general  nonlinear  constrained 
optimization  problem  can  be  stated  mathematically  as: 


Minimize: 

F{X) 

objective  function 

(2.91) 

Subject  to: 

9j{^)  <  0 

j  =  l,rn 

inequality  constraints 

(2.92) 

hk{X)  =  0 

k  =  1,1 

equality  constraints 

(2.93) 
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xj  <  Xi  <  x^ 


i  —  1 ,  n 


side  constraints 


(2.94) 


where 


The  design  variables  are  defined  as  Xi  to  Xn. 


(2.95) 


Most  optimization  algorithms  start  from  an  initial  guess  and  are  iterated  upon  to  find 
the  optimal  solution.  The  most  common  form  of  the  iterative  procedure  is 


X<}  =  +  0*8“  (2.96) 

where  q  is  the  iteration  number,  S  is  the  search  direction  in  the  design  space  and  0*  is  the 
optimal  step  length  in  that  direction  . 

2.2.1 .1  Unconstrained  Optimization.  From  calculus,  we  know  that  for  F(X) 
to  be  a  minimum,  the  gradient  of  F(X)  must  be  equal  to  zero.  This  condition  is  necessary 
to  say  the  point  is  a  local  minimum,  but  does  not  guarantee  a  global  minimum.  Also  from 
calculus,  we  know  that  if  the  second  derivative  of  the  function  with  respect  to  all  variables  is 
positive  definite,  that  point  is  a  minimum.  This  leads  to  the  definition  of  the  Hessian  matrix. 
The  Hessian  matrix  is  a  matrix  of  second  partial  derivatives  of  the  function  with  respect  to 
the  design  variables.  See  Equation  (2.97).  If  the  Hessian  is  positive  definite  and  the  gradient 
is  equal  to  zero,  this  ensures  a  local  minimum.  The  design  is  a  global  minimum  only  if  the 
Hessian  is  positive  definite  for  all  possible  values  of  the  design  variables  within  the  design 
space. 
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r  a^F(X) 

ax'{ 

d^F(X) 

dXiax2 

a'^Fix)  1 
axiaxn 

H  = 

a^F{X) 
9.^2  9^1 

... 

^xf 

a'^Fix) 

aX2dXn 

(2.97) 

a^F(X) 

a'^F(X) 

a'^F(x) 

L  aXndXi  dX„dX2  dXl  -I 

2.2.12  Constrained  Optimization.  For  constrained  optimization,  it  isn’t 
necessary  that  the  gradient  of  the  objective  vanish  at  the  optimum.  This  type  of  optimization 
problem  works  to  reduce  the  objective  without  violating  the  constraints.  This  is  done  by 
finding  a  search  direction.  A  direction  that  will  reduce  the  objective  is  called  usable.  A 
feasible  direction  is  defined  as  one  that  wiU  not  violate  a  constraint  if  a  small  step  is  taken.  An 
active  constraint  is  defined  as  one  that  is  within  some  tolerance  of  zero.  A  search  direction 
is  found  that  is  feasible  with  respect  to  all  active  constraints.  As  stated,  the  search  direction 
must  be  both  usable  and  feasible. 


Usable  Direction;  VF(A’)  •  5  <  0  (2.98) 

Feasible  Direction;  \^gj{X)  ■  S  <0  all  j  for  which  gj{X)  =  0  (2.99) 

The  optimum  of  a  constrained  optimization  problem  can  be  defined  as  meeting  three 
criteria  called  the  Kuhn-Tucker  Necessary  Conditions.  These  conditions  are; 


1.  X*  is  feasible 

(2.100) 

2.  A,-^,(X*)=0 

j  =  l,m 

A,>0 

i 

Y,  Afc+„,  V/ifc(X*)  =  0 

(2.101) 

3.  VF{X*)  + 

A,  V5,(X-)  + 

(2.102) 

k=l 

Xj  >  0  unrestricted 
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where  X*  is  the  vector  of  optimum  design  variables,  g  is  the  vector  of  inequality  constraints, 
F  is  the  objective,  h  is  the  vector  of  equality  constraints  and  the  A’s  are  Lagrange  multipliers. 
The  first  Kuhn-Tucker  condition  (2.100)  states  that  the  optimal  solution  must  satisfy  all  the 
constraints.  The  second  Kuhn-Tucker  condition  (2.101)  says  that  if  the  constraint  is  not  equal 
to  zero,  then  the  Lagrange  multiplier  must  be  equal  to  zero.  The  last  Kuhn-Tucker  condition 
states  that  the  objective  gradient  at  optimum  must  be  in  the  subspace  spanned  by  the  active 
constraint  gradients.  Also,  for  the  Kuhn-Tucker  conditions  to  be  necessary  for  a  constrained 
optimum,  the  X*  must  be  a  regular  point,  (i.e.,  a  design  where  the  active  constraint  gradients 
are  linearly  independent).  If  the  objective  function  and  all  constraint  surfaces  are  convex, 
then  the  design  space  is  said  to  be  convex  and  the  necessary  Kuhn-Tucker  conditions  are  also 
sufficient  to  guarantee  that  the  optimum  reached  is  a  global  optimum. 

2 .2 .2  History  of  the  Mixed  Numerical  Method.  This  section  addresses  the  numerical 
approach  to  designing  a  mixed  H2IH00  controller.  This  method  approaches  an  optimal  fixed 
order  solution.  The  controller  from  this  design  method  wiU  not  satisfy  the  constraint  exactly, 
but  is  as  close  as  numerically  possible  within  some  tolerance. 

Ridgely  [Rid91]  developed  a  numerical  solution  using  a  new  performance  index.  This 
was  based  on  a  connection  between  the  central  Hoo  controller  [DGKF89]  and  the  mixed 
controller.  This  performance  index  is 

=  (1  -  ^)||7;,„||^  -f  g  tr[Q^CjC,]  (2.103) 

The  value  p,  was  varied  from  one  to  zero,  resulting  in  a  mixed  controller.  When  p 
equals  1,  the  central  Hoo  controller  is  the  solution.  The  DFP  method  was  used  in  optimizing 
this  set-up.  There  were  several  drawbacks  to  this  method: 

1.  Amount  of  computational  time  required  was  on  the  order  of  one  week  for  one  point  on 
the  mixed  H2IH00  curve,  since  each  p  design  was  optimized  to  a  point  on  the  mixed 
curve. 
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2.  The  order  of  the  controller  achieved  by  this  method  is  the  order  of  the  mixed  plant 
together  with  all  of  its  weights. 

3.  This  method  does  not  allow  for  singular  Hoo  constraints. 


Walker  developed  a  numerical  approach  which  restated  the  mixed  H2 /Hoo  optimization 
problem,  where  the  two-norm  is  the  objective  and  the  infinity-norm  is  the  constraint.  The 
initial  guess  for  the  optimization  process  was  the  optimal  controller,  since  it  was  a  relatively 
easy  point  to  calculate.  Then  steps  in  the  infinity-norm  constraint  were  made  to  develop  the 
a*  versus  7*  curve.  Inequality  constraints  were  used  to  develop  this  performance  index: 


min 

K  stabilizing 


subject  to 


||re<i(X)||oo-7<0 


This  method  saved  time  compared  to  Ridgely’s  method  since  each  controller  found  was 
a  point  on  the  mixed  curve.  Walker  initially  used  the  DFP  method  to  minimize  the  performance 
index.  To  prevent  non-stabilizing  controllers,  the  performance  index  was  set  to  an  artificially 
high  value  if  an  unstable  closed-loop  controller  was  calculated.  Thus,  that  controller  would 
not  be  accepted.  Initially,  equality  constraints  were  addressed,  but  due  to  numerical  drawbacks 
Walker  then  turned  to  using  an  inequality  constraint  approach.  The  final  algorithm  used  to 
minimize  the  performance  index  was  the  MATLAB  version  of  the  SQP  technique. 

2.2.5  SQP  Optimization.  The  Sequential  Quadratic  Programming  method  uses  a 
technique  where  the  search  direction  is  found  by  solving  a  quadratic  objective  function  and 
linear  constraints.  The  search  direction  is  found  by  creating  a  quadratic  approximation  to  the 
objective  function.  The  general  problem  set-up  is  as  follows: 

Minimize:  g(5)  =  F(X) -h  [VF(A:)]^5 -f  (2.104) 
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Subject  to:  Vgj{X)*S  +  8jgj{X)  <0  j  =  l,m  (2.105) 

Xhk{X)*S  +  Shk{X)  =  0  k  =  l,l  (2.106) 

The  design  variable  for  this  subproblem  is  S,  the  search  direction.  The  matrix  B  is  updated 
through  iterations  to  approach  the  Hessian  matrix.  The  ^’s  ensure  that  the  linearized  constraints 
do  not  cut  off  the  feasible  region  of  the  design  space.  These  are  problem  dependent  and  should 
be  chosen  as  close  to  1  as  possible. 


8j  =  l  ifgj{X)<0  (2.107) 

8j=8  ifgj{X)>0  (2.108) 

0  <  ^  <  1  (2.109) 

Once  the  search  direction  is  found,  a  one  dimensional  search  is  performed  with  an  exterior 
penalty  function  added. 


m  I 


$  =  F{X)  -f-  Uj  {max  [0,5rj(A:)]}  +  Y,  ^m+k  [hk{X)] 

j=i  k=i 

(2.110) 

X  =  X^-^  -f  ps 

(2.111) 

=  l^il  i  =  l,m  +  / 

first  iteration  (2.112) 

Uj  =max[|Aj|,i  (u' +  |Aj|) 

next  iterations  (2.113) 

The  variable  u  is  a  vector  penalty  parameter.  The  way  this  problem  is  set  up,  the  one 
dimensional  search  is  well  conditioned  and  the  step  of  ^  =  1  is  a  good  step  size  to  make 
once  the  matrix  B  approaches  the  Hessian  matrix.  Now  that  the  first  step  has  been  taken,  the 
B  matrix  needs  to  be  updated.  The  B  matrix  improves  the  quadratic  approximation  for  the 
search  direction  subproblem.  There  are  several  methods  used  to  update  B.  The  method  used 
by  MATLAB  is  the  Broydon-Fletcher-Shanno-Goldfarb  (BFGS)  update  formula  [Van84]. 


B* 


^  Bpp^B  ^  rjT)'^ 
p^Bp  p  •  rj 


(2.114) 
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p  = 

rj  =  Oy  +  {1-  6)Bp 
y  =  A")  -  A'' 

m  / 

«Jf,A)  =  F{X)-\-^Xigi{X)+-£X„ 

j=l  k=l 

l.Oif  p  •  ?/  >  0.2p^Bp 
<  O-2/^P 

This  update  formula  maintains  positive  definiteness  and  is  called  a  variable  metric 
method.  This  update  method  also  allows  for  infeasible  designs  with  no  other  additional 
methods  required.  Since  the  solution  is  allowed  to  approach  from  the  infeasible  design  space, 
the  resulting  controller  may  be  non-stabilizing.  The  previous  performance  index  and  problem 
set-up  discussed  in  Section  2.2.2  can  be  used,  but  an  additional  stability  constraint  must  be 
added. 

2.2.4  Two-Norm  Gradients.  The  gradient  of  the  two-norm  at  H2  optimal  was 
defined  by  Ridgely  [Rid91].  He  used  the  definition  of  the  two-norm 


(2.115) 

(2.116) 

(2.117) 

+khkiX) 

(2.118) 

] 

(2.119) 

\\T.w\\l  =  tr  [Q2CJCJ  (2.120) 

The  H2  problem  can  be  set  up  as  an  optimization  problem  by  forming  a  Lagrangian.  This 
Lagrangian  has  the  objective  function  as  the  definition  of  the  two-norm.  The  constraint  in  the 
Lagrangian  is  the  Lyapunov  equation  that  must  be  solved  to  find  a  real,  positive  semidefinite 
solution  for  Q2  in  Equation  (2.120).  The  Lagrangian  is  then 


C  —  tr  Q2CJCz  -f-  tr  {A2Q2  +  Q2A2  -f 


(2.121) 


Taking  the  partial  derivatives  of  this  equation  with  respect  to  all  the  unknown  variables 


yields 


dC, 

dAc 

dC 

dBc 

dC 

dC 

W2 

dc 


—  ^[^uQu  +  ^2Q7\ 


=  2[XlQ,Cl  +  X^Ql^Cl  + 


—  2\,B^^X\Q\2  +  B^^Xx^Qi  +  RJ2Q12  +  R2CCQ2] 
=  A'^  X  +  XA2  + 

=  A2Q2  +  <52-^2  +  ^wR'w 


(2.122) 

(2.123) 

(2.124) 

(2.125) 

(2.126) 


The  process  in  which  to  calculate  the  H2  gradient  with  respect  to  a  specific  controller  is  as 
follows 


1.  Set  (2.125)  and  (2.126)  equal  to  zero  and  solve  for  X  and  Q2 

2.  Compute  and  ^  from  ( 2.122),  ( 2.123),  and  ( 2.124),  respectively 

3.  Rearrange  the  gradient  by 


dC 


(2.127) 


A  method  to  solve  the  derivative  of  the  two-norm  at  a  sub-optimal  point  was  derived 
by  Canfield  [Can94].  He  showed  that  the  derivative  of  the  two-norm  can  be  derived  through 
a.  Direct  Method.  The  derivation  begins  with  the  definition  of  the  two-norm. 


\\T.4\l  =  tr[Q2CjC,] 


where  (^2  is  the  solution  to  the  following  Lypaunov  equation 


^2^2  +  <?2^2  +  RwB^  —  0 


(2.128) 


(2.129) 
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The  closed-loop  matrices  A2,  and  are  explicit  functions  of  the  controller  matrices 
Ac,  Be,  and  Cc  which  are  parameterized  by  the  design  vector  P.  Equation  (2.128)  can  be 
differentiated  with  respect  to  the  design  vector  P. 

\\T,X  =  tr[Q2'CjC,  +  QCJ'C,  -f  Q2CJCJ]  (2.130) 


Since  Q2  is  an  implicit  function  of  P,  the  Lyapunov  equation  (2.129)  can  be  differentiated 
with  respect  to  each  design  variable.  This  is  given  by 

^2Q2  +  A2Q2  +  Q2  +  Q2^  +  BjB^  •+•  B^B^  =  0  (2.131) 


This  can  then  be  rearranged  to  form  another  Lyapunov  equation: 


A2Q2'  +  Q2A2  +  {A2'Q2  +  Q2^  +  BJB^  -f  By^B^  )  =  0  (2.132) 


Then  (2.132)  can  be  solved  for  Q2  for  each  design  variable.  The  derivative  of  the  two-norm 
is  then  found  by  substituting  Q2  into  (2.130).  This  method  was  not  used  to  calculate  the  H2 
gradients  at  sub-optimal  mixed  H2/H00  designs  since  a  Lyapunov  equation  would  have  to  be 
solved  for  every  design  variable  for  every  gradient.  This  would  be  too  expensive.  Rather, 
Canfield  also  derived  an  Adjoint  Method  for  the  gradient  of  the  two-norm. 

The  Adjoint  Method  begins  with  the  definition  of  the  two-norm  which  will  be  called  ^ 
for  this  derivation. 

^{p,  Q2)  =  tr[Q2Cj{p)C,{p)]  (2.133) 

Take  the  first  variation  of  (2.133),  treating  p  and  Q2  as  independent  variables. 


«(?,  <?.)  =  E  +  E  E 


•jk 


(2.134) 


However,  variations  6Q2jk  are  not  independent  of  Spu  rather,  they  are  related  through  the 
Lyapunov  equation  (2.129).  An  identity  can  be  formed  by  post-multiplying  (2.129)  by  the 
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adjoint  matrix  X  and  taking  the  trace. 


Qii^)  —  ^^[(-42^2  +  Q2^  +  ^•u}B^)X]  —  0  (2.135) 


The  equation  (2.135)  will  be  called  the  adjoint  equation.  The  variation  of  the  adjoint  equation 
is  then 


aT(p,  Q2.X)  _  ^  5T(p,  X) 

- - hi  +  1^1^ - =  0  (2.136) 


dQi 


j  k 

where  X  is  not  free  to  vary,  but  instead  is  defined  as  the  matrix  that  satisfies 


drip,Q2,x)  _  d^ip,Q2) 

dQ2  dQ2 

which  when  combined  with  equation  (2.133)  and  (2.135)  can  be  written  as 


(2.137) 


A^X-^  XA2  +  C,Cj  =  0 


(2.138) 


Therefore  (2.138)  is  a  Lyapunov  equation  for  X  that  does  not  depend  on  any  derivatives 
with  respect  to  p  (i.e.,  it  is  only  solved  once).  Substituting  equation  (2.137)  into  (2.136)  the 
following  is  found 


EE 


d<l(p,Q2,X) 

8Q2„ 


=  E 

i 


dT{p,Q2,X) 

dpi 


6pi 


(2.139) 


which  provides  the  necessary  relationship  to  determine  how  small  variations  Sp  of  the  design 
parameters  cause  corresponding  small  variations  6Q2  of  the  Lyapunov  matrix  Q2.  Next  the 
total  variation  can  be  expressed  by  substituting  equation  (2.139)  into  (2.134). 


i  dpi  dpi 


(2.140) 


2-28 


where  the  terms  in  brackets  are  recognized  as  the  design  derivatives  that  we  seek 


a||r..||2  mp.Q-^)  arM2,x) 

dp  dp  dp 


(2.141) 


where 


dcj 


dpi 


dT{p,Q,,X) 

dpi 


,dA 


dpi 

dAi:  dB,. 


dpi 


— Q2  +  Q2-^ - 1-  — B^  +  B, 


dp, 


dpi  dpi 


dBl 

’  dpi 


)^] 


(2.142) 

(2.143) 


and  where  Q2  is  found  by  solving  the  Lyapunov  equation  (2.129)  and  X  is  found  by  solving 
the  Lyapunov  equation  (2.138). 

Equation  (2.143)  represents  the  change  in  the  two-norm  with  respect  to  the  design 
variables  in  vector  form.  Equation  (2.143)  can  be  used  to  show  the  change  in  the  two-norm 
with  respect  to  the  design  variables  in  matrix  format  (i.e.,  Ac,  Be  and  Cc). 


M(Ac, 

Be, 

Cc,Q2) 

dAe 

M(Ae, 

Be, 

Cc,Q2) 

dBe 

d^{Ac, 

Be, 

Cc,Q2) 

dCe 

dT^AcBc 

,Ce, 

Q2,X) 

dAc 

dT{Ac,Bc 

,Ce, 

Q2,X) 

dBc 


dT{Ac,Bc,Cc,Q2,X) 

dCc 


{Q2 

Q2CI 


0  D. 


0  D. 


Cz  + 

) 


(2.144) 

(2.145) 

(2.146) 


0  0 
0  1 

0  0 
Cy2  0 
0 

Byyj 


Q2  +  Q2 

Q2  +  Q2 


0  0 
0  1 


0  Cy2 
0  0 


+ 


-|-  B,ui 


0  D 


yw 


0  Bu2 
0  0 


Q2  +  Q2 


0  0 

Bu2  0 


(2.147) 

(2.148) 


(2.149) 
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By  simplifying  Equations  (2.144)  through  (2.149)  and  substituting  them  into  (2.142)  through 
(2.143)  and  then  substituting  the  result  into  (2.141),  in  turn,  the  final  equations  for  the  derivative 
of  the  two-norm  for  a  suboptimal  mixed  H2fHoo  design  are  found.  The  result  can  also  be 

I  as 

=  2[X^(5i2  +  -^2^2]  (2.150) 

=  21a:5(J,cJ  +  (2.151) 

“  2[B^XiQi2  +  +  ^12^12  (2.152) 

Therefore,  the  right  side  of  (2.150)  through  (2.152)  is  the  right  side  of  (2.122)  through  (2.124). 

The  two-norm  can  be  defined  for  both  a  stable  and  unstable  closed-loop  system.  The 
two-norm  is  defined  for  an  unstable  transfer  function  as  long  as  the  system  has  no  poles  on 
the  imaginary  axis  and  is  strictly  proper.  If  the  system  is  stable  the  two-norm  is  calculated 
by  the  methods  discussed  in  Section  2.1.1.  If  the  system  is  unstable,  it  can  be  broken  into  its 
stable  and  anti-stable  parts.  The  two-norm  of  an  anti-stable  transfer  function  is  the  two-norm 
of  its  complex  conjugate  transpose,  which  is  stable.  Therefore,  the  two-norm  of  an  antistable 
transfer  function  can  be  found  by  setting  the  state  space  matrices  of  the  antistable  system  equal 
to  A  =  B  =  and  C  =  —B^.  Then  the  two-norm  can  be  calculated  by  the  methods 
discussed  in  Section  2.1.1.  Finally,  through  orthogonality,  ||G||2  =  ||nG||2-l- ||n-‘-G||2  where 
IIG  is  the  stable  projection  and  is  the  antistable  projection.  The  two-norm  and  therefore 
the  gradients  must  be  calculated  for  stable  and  unstable  closed-loop  systems. 


expressed  in  matrix  forr 


d\\Z 


zwl\2 


dAr 


m 


2 

zwll2 


dBr 


m 


zw\\2 


dCr 


2.2.5  Infinity-Norm  Gradients.  The  infinity-norm  gradient  can  be  calculated  by 
an  analytic  method.  This  method  is  called  a  singular  value  sensitivity  approach,  which  was 
developed  by  Geisy  and  Lim  [GL93].  Assuming  the  maximum  singular  value  of  Ted  at  X  has 
a  single  peak,  the  derivatives  of  the  infinity  norm  can  be  evaluated  as: 


^||Ted[|oo 

dxi 


=  3?e 


/aTerf(a;o)\ 

V  dxi  ) 


(2.153) 
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where  ui  and  vi  are  the  singular  vectors  associated  with  the  maximum  singular  value  of  Ted. 
The  Uo  is  the  frequency  at  which  the  peak  occurs,  Xq  is  the  current  guess  of  X  and  3?e  denotes 
the  real  part.  The  derivative  of  the  transfer  function  with  respect  to  the  design  variables  can 
be  written  as 

dTed{i^o)  ^  (i^o  ~  *^00)  +  ^ed] 

dxi  V 

Aq 

There  are  several  advantages  to  this  method  over  finite  differencing: 

1.  The  analytical  method  is  quicker  since  there  is  only  one  gradient  calculation  per  design 
variable.  The  finite  difference  method  requires  an  infinity-norm  calculation  with  the 
design  variables  perturbed  forward  and  backward  and  then  an  average  is  taken  between 
the  two.  The  finite  difference  method  therefore  requires  twice  as  many  calculations  as 
the  analytical  derivative. 

2.  The  accuracy  of  the  analytic  gradient  is  better  since  in  finite  difference  the  solution 
is  based  on  differences  between  very  small  numbers  which  can  result  in  errors  due  to 
round  off. 

3.  The  Hcc  norm  can  be  identified  as  a  piecewise  continuous  function  since  the  norm 
definition  is  the  peak  of  the  maximum  singular  value  plot.  The  finite  difference  ap¬ 
proach  perturbs  the  controller  and  then  finds  the  norm.  The  frequency  where  the  peak 
is  found  could  change  and  the  gradient  at  a  different  frequency  is  usually  not  the  same 
as  the  previous  peak.  Therefore,  a  central  difference  method  could  result  in  an  incorrect 
gradient.  The  analytic  method  does  not  have  this  problem,  since  it  provides  information 
only  on  the  unperturbed  transfer  function. 

However,  the  sensitivity  method  has  some  disadvantages: 


(2.154) 
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1.  The  method  has  difficulty  converging  when  the  peaks  are  located  at  more  than  one 
frequency.  This  results  in  a  discontinuous  gradient.  This  is  a  problem  for  both  finite 
difference  and  the  analytic  gradient  method. 

2.  The  method  depends  on  a  search  range  over  frequency.  Previous  knowledge  of  the 
shape  of  the  singular  value  plot  can  be  used  to  specify  this  range.  A  check  needs  to  be 
made  that  the  range  given  includes  the  maximum  singular  value. 

3.  The  fineness  of  the  frequency  search  can  dictate  the  accmacy  of  the  Hoo  norm. 
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Ill.  Optimization  Code  Improvements 

This  section  is  divided  into  two  parts.  Section  3.1  discusses  improvements  made  to 
the  code.  For  each  improvement  the  problem  definition  and  solution  are  given.  Section  3.2 
discusses  an  implementation  of  the  improvements  for  a  simple  SISO  example. 

3.1  Code  Improvements 

There  were  eight  general  improvements  made  to  the  mixed  H2  / Hoo  optimization  code. 
They  are  listed  in  order  of  importance  relative  to  the  improvement  in  run  time  and  accuracy 
of  the  optimization  code. 

3 .1 .1  Peak  Identification. 

Problem  Definition:  The  infinity-norm  of  Ted  is  defined  as  the  supremum  over  all  frequency 
of  the  maximum  singular  value  of  Ted.  There  can  be  several  peaks  in  the  singular  value 
plot  of  a  system.  It  was  found  that  the  optimization  routine  could  not  handle  the  case  when 
two  or  more  peaks  had  approximately  the  same  singular  value.  The  routine  would  constrain 
the  peak  at  one  frequency  while  the  other  frequency  was  not  constrained.  This  caused  the 
optimizer  to  switch  back  and  forth  between  the  peak  values  and  not  be  able  to  complete  the 
optimization  iteration  where  the  multiple  peaks  occurred.  This  is  an  identified  limitation  of 
using  the  analytic  gradient  method  discussed  by  Geisy  and  Lim  [GL93]. 

Problem  Solution:  The  solution  to  this  problem  was  to  identify  all  peaks  (and  their  frequencies) 
which  were  within  a  specified  tolerance  of  the  maximum.  The  peaks  were  identified  by 
evaluating  the  change  in  slope  of  the  singular  value  plot.  The  peaks  were  then  separately 
constrained.  The  singular  value  at  a  frequency  of  zero  rad/sec  was  automatically  constrained 
and  a  frequency  search  range  of  10“®  to  10^  rad/sec  was  used.  A  frequency  band  was  placed 
around  each  peak  to  allow  for  movement  of  the  peak  as  the  controller  design  changed  in  the 
optimization  routine.  A  function  was  implemented  that  allowed  for  a  ±  0.3  rad/sec  band  about 
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Singular  Value  Plot 


Figure  3 . 1  Singular  Value  Plot  (7=20) 


the  peak  and  also  prevented  the  bands  from  overlapping.  Two  examples  of  singular  value 
plots  identifying  peaks  and  bands  are  shown  in  Figures  3.1  and  3.2. 

The  frequency  search  for  peak  identification  was  done  before  and  after  each  7  step. 
The  search  before  the  optimization  was  done  to  set  up  the  bands  once  a  peak  was  identified. 
The  search  after  the  optimization  was  to  check  that  another  peak  had  not  risen  above  the 
constrained  7  level  outside  the  bands  previously  set.  If  a  peak  was  identified  outside  the  bands 
of  the  optimization,  the  bands  from  before  the  optimization  and  after  the  optimization  would 
be  combined.  The  bands  were  sorted  and  the  frequencies  were  compared  in  log  scale.  If 
the  difference  between  each  band  was  less  than  0.05  in  log  scale,  then  that  band  would  be 
removed.  This  was  a  way  to  prevent  overlapping  and  duplication  of  bands.  Once  the  bands 
were  reordered,  the  optimization  sequence  would  begin  again  with  the  new  bands  and  the  last 
design  vector  found.  The  optimization  sequence  would  not  end  until  the  H^o  constraint  was 
completely  satisfied. 

3.1.2  Gradient  Calculation.  The  following  five  subsections  wiU  discuss  the  cor¬ 
rections  made  to  the  gradient  calculations.  They  are  listed  in  order  of  importance.  First,  the 
correction  to  the  Hoo  gradient  is  presented  and  then  the  correction  made  to  the  H2  gradient 
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Singular  Value  Plot 


Figure  3 .2  Singular  Value  Plot  (7= 1 .3  6) 


follows.  Finally,  introduction  of  an  analytic  means  to  take  the  derivative  of  the  stability 
constraint  (previously  done  using  finite  differencing)  is  given. 

3. 1.2.1  Hoo  Gradient  Modification. 

Problem  Definition:  The  analytic  method  for  calculating  the  Hoo  gradient  is  faster,  more 
accurate  and  better  conditioned  for  finding  a  numeric  answer  than  the  finite  difference  method. 
However,  when  comparing  the  results  of  these  two  gradient  methods,  the  gradients  with 
respect  to  each  design  variable  should  be  close  in  value.  When  comparing  the  Hoo  gradients, 
the  results  were  significantly  different.  The  gradients  were  opposite  in  sign  in  some  cases. 
Another  indication  of  a  gradient  problem  was  that  the  SQP  optimization  routine  was  regularly 
giving  negative  step  sizes.  This  indicated  the  search  direction  the  optimizer  was  calculating 
might  be  in  error.  This  led  to  an  investigation  of  the  analytic  gradient  calculation  for  the 
infinity-norm. 

Problem  Solution:  The  background  for  the  infinity-norm  calculation  is  presented  in  Sec¬ 
tion  2.2.5.  The  two  equations  that  define  the  analytic  Hoo  gradient  are 
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^||rerf||oo 

dxi 


H  ( dTedi'^o)^ 

V  dx,  ) 

V\ 

(3.1) 


and 


dTedij^o) 

d 

Ce  (i^^O  ■“  .^oo  )  ^  Bd-\-  T>ed 

dxi 

X, 

dXi 

The  important  part  of  (3.2)  is  that  it  defines  the  derivative  of  the  closed  loop  transfer  function 
Ted  with  respect  to  each  value  in  the  controller  state  space  (i.e.,  the  design  variables  for  a 
numerical  optimization  solution).  Using  these  two  equations,  the  derivative  of  the  infinity- 
norm  of  the  closed  loop  Ted  transfer  function  can  be  found  with  respect  to  each  design  variable. 

The  original  set-up  of  the  mixed  problem  was  to  put  the  controller  state  space  into 
modal  form.  Modal  form  places  the  real  parts  of  the  eigenvalues  of  the  controller  A  matrix 
on  the  diagonal,  with  the  imaginary  parts  of  the  eigenvalues  on  the  superdiagonal  and  the 
negative  of  the  imaginary  parts  of  the  eigenvalues  on  the  subdiagonal.  An  example  of  a  modal 
transformation  is 


A 


14  5  7 
3  8  2  9 
3  5  7  1 
0  3  6  8 


Amodal 


18.10 

0.000 

0.000 

0.000 

0.000 

0.117 

0.000 

0.000 

0.000 

0.000 

2.886 

3.662 

0.000 

0.000 

-3.662 

2.886 

(3.3) 
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This  modal  form  significantly  reduces  the  amount  of  design  variables  needed  to  be 
found  in  an  optimization  routine.  The  higher  the  order  of  the  controller,  the  more  the  savings 
in  design  variables.  The  reduction  in  number  of  design  variables  with  the  modal  form  comes 
from  the  fact  that  only  the  diagonal  terms  and  superdiagonal  terms  in  the  A  matrix  need  to 
be  found.  The  subdiagonal  terms  do  not  need  to  be  calculated  since  they  are  just  the  negate 
of  the  superdiagonal.  The  state  space  B  and  C  matrices  for  the  controller  are  assumed  to  be 
fully  populated.  An  example  of  a  fourth  order  controller  A  matrix  is 


Amodal 


Xi 

X5 

0 

0 

-X5 

X2 

Xe 

0 

0 

-Xe 

X3 

X7 

0 

0 

-X7 

X4 

(3.4) 


Notice  for  this  example  that  putting  the  A  matrix  into  modal  form  reduces  the  required  number 
of  design  variables  in  A  from  16  down  to  7. 

The  derivative  as  defined  in  (3.2)  must  be  taken  with  respect  to  every  design  variable. 
Previously,  the  gradients  with  respect  to  the  subdiagonal  terms  were  not  factored  into  the 
overall  gradient.  This  caused  the  gradients  to  be  in  error.  This  error  caused  examples  of 
the  search  direction  generated  by  the  SQP  optimization  routine  to  be  in  opposite  directions 
from  what  the  finite  difference  gradient  calculated.  Adding  the  gradient  with  respect  to  the 
subdiagonal  terms  of  the  A  matrix  corrected  this  problem,  and  the  finite  difference  matched 
very  well  with  the  analytic  gradients  for  the  infinity-norm. 


3.1 .2.2  H2  Modification  for  Modal  Form. 

Problem  Definition:  The  background  of  the  H2  gradient  is  given  in  Section  2.2.4.  This 
development  showed  that  the  resulting  gradients  are  a  function  of  taking  the  partial  derivative 
of  the  Lagrangian  function  with  respect  to  each  of  the  controller  state  space  matrices.  Since 
there  is  a  partial  with  respect  to  the  controller  A  matrix  in  the  definition  of  the  H2  gradient. 
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the  subdiagonal  gradient  terms  were  neglected,  similar  to  the  Hoo  gradient  problem.  The 
following  shows  the  Lagrangian  equation  and  the  partial  of  the  Lagrangian  with  respect  to  Ac, 
(i.e.,  the  controller  A  matrix): 


C  =  tr  [gsCjc, 


+  tr  {A2Q2  +  Q2A2  + 


(3.5) 


—  =  2iX^,Q,2  +  X2Q2)  (3.6) 

The  problem  with  the  H2  gradients  was  identified  in  the  same  way  as  the  Hoo  gradients.  The 
H2  gradients  were  significantly  in  error  from  finite  difference  without  the  subdiagonal  terms. 


Problem  Solution:  Adding  the  derivative  with  respect  to  the  subdiagonal  terms  to  the  equation 
for  the  partial  of  the  Lagrangian  with  respect  to  Ac  corrected  this  problem. 

3.1 .23  H2  Modification  for  Sub-Optimal  Design. 

Problem  Definition:  Comparison  of  the  finite  difference  to  the  analytical  two-norm  gradients 
for  each  design  variable  showed  a  difference  of  a  factor  of  two.  This  lead  to  a  further  investi¬ 
gation  of  the  analytic  two-norm  gradient.  Through  calculation  of  the  necessary  conditions  it 
was  found  that  there  was  a  factor  of  two  in  each  of  these  derivatives.  This  factor  was  ignored 
in  the  two-norm  gradient  equations,  since  this  calculation  was  considered  at  the  H2  optimal 
point  and  therefore  the  necessary  conditions  would  be  equal  to  zero. 


Problem  Solution:  If  the  mixed  design  is  not  at  the  H2  optimal  point,  the  mixed  design  is  by 
definition  sub-optimal.  Therefore,  the  gradients  of  the  Lagrangian  of  H2,  when  completing  a 
mixed  H2/ Hoo  optimization  design,  should  contain  the  factor  of  two  for  each  derivative  with 
respect  to  the  controller  state  space  matrices.  The  two-norm  derivative  at  sub-optimal  designs 
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can  be  written  as 


OL 

dA, 

=  2  *  [X^Qi2  +  X2Q‘^ 

(3.7) 

dC 

dB, 

=  2  *  [XlQ.Cl  -b  X2QI2CI  +  X^^V,2] 

(3.8) 

dC 

dC, 

—  2  +  B^^X\2Q2  +  Ri2Qi‘2  +  R2CcQ2\ 

(3.9) 

Solve  for  X  : 

[A'^X  -|-  XA2  -h  —  0 

(3.10) 

Solve  for  Q2  ■ 

[•^2Q2  +  Q2-^2  +  ^w^w]  —  0 

(3.11) 

3. 1.2. 4  H2  Gradient  Modification  if  Unstable. 

Problem  Definition:  The  SQP  optimization  method  is  allowed  to  begin  with  an  initial  guess 
that  is  in  the  infeasible  region.  SQP  also  allows  a  feasible  solution  to  go  into  the  infeasible 
region  during  the  optimization  process.  A  non-stabilizing  controller  is  defined  when  the 
stability  constraint  is  violated,  and  therefore  the  design  is  in  the  infeasible  region.  Since  it  is 
possible  that  the  controller  could  be  non-stabilizing,  the  code  should  be  able  to  calculate  the 
two-norm  for  this  case.  The  previous  code  did  a  stable/anti-stable  projection  to  calculate  the 
two-norm,  but  did  not  consider  this  projection  for  the  gradient  of  the  two-norm. 

Problem  Solution:  An  addition  was  made  to  the  code  to  calculate  the  H2  gradient  when 
the  controller  was  identified  as  non-stabilizing.  This  gradient  was  calculated  using  finite 
differencing.  The  section  of  the  code  that  calculates  the  gradient  using  finite  differencing  is 
called  whenever  the  real  part  of  the  maximum  eigenvalue  of  the  closed  loop  A2  matrix  of  the 
H2  problem  is  positive. 

3. 1.2. 5  Stability  Gradient  Modification. 

Problem  Definition:  The  previous  method  calculated  the  stability  gradient  using  the  finite 
difference  method.  This  method  became  costly  due  to  the  fact  that  for  every  design  variable, 
two  eigenvalue  evaluations  must  be  done.  This  must  be  done  for  every  derivative  calculation  in 
the  gradient.  A  method  was  derived  to  find  the  gradients  of  the  stability  constraint  analytically. 
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Problem  Solution:  The  stability  constraint  is  defined  as  requiring  the  real  part  of  the  maximum 
eigenvalue  of  the  closed-loop  A2  matrix  to  be  negative.  The  derivative  of  this  constraint  is  the 
derivative  of  the  eigenvalue.  The  eigenvalue  of  the  closed-loop  A2  matrix  can  be  expressed 
as 

A  =  u^A2V2  (3.12) 

where  U2  and  V2  are  the  left  and  right  eigenvectors,  respectively,  associated  with  the  eigenvalue 
A.  The  stability  constraint  is 

max3?e(A)  <  0  (3.13) 

so  the  stability  constraint  gradient  is 

=  3Je[uf  ^^^2]  (3.14) 

The  mixed  H2IH00  A2  matrix  is  an  explicit  function  of  the  design  variables.  Taking  the 
derivative  of  A2  with  respect  to  each  design  variable,  and  then  pre-  and  post-  multiplying  it 
by  the  normalized  eigenvectors  produced  the  analytic  stability  gradient.  This  allowed  one 
eigenvalue  evaluation  per  stability  gradient  instead  of  2n+l  evaluations,  where  n=number  of 
design  variables.  Although  eigenvalue  calculations  are  not  expensive,  the  cost  increases  with 
an  increase  in  the  number  of  design  variables.  Note  that  this  development  assumes  that  the 
maximum  eigenvalue  is  not  repeated.  If  it  were  a  repeated  eigenvalue,  further  modifications 
would  need  to  be  made.  These  were  not  considered  here. 

3.13  Penalty  Function. 

Problem  Definition:  SQP  allows  infeasible  designs  in  the  optimization  process;  therefore, 
a  design  (i.e.,  controller)  may  be  non-stabilizing.  If  the  controller  is  non-stabilizing,  the 
calculation  of  the  two-norm  gradient  may  be  very  expensive.  The  expense  is  due  to  the  finite 
differencing  routine  that  was  added  to  calculate  the  two-norm  gradient  when  the  stability 
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constraint  was  violated.  Therefore,  it  would  be  advantageous  to  prevent  non-stabilizing 
controller  designs  while  in  the  optimization  process. 


Problem  Solution:  The  solution  to  this  problem  was  to  add  a  penalty  to  the  objective  function 
which  was  a  function  of  the  stability  constraint.  The  penalty  function  (j)  added  is 

=  [A3  *  max(real(eig(v42)))]^  (3.15) 

where  A3  is  the  scale  factor  used  to  prevent  unstable  designs.  This  will  be  discussed  in  the 
scaling  section. 

If  the  controller  is  non-stabilizing,  the  gradient  of  the  two-norm  must  also  have  the 
penalty  function  gradient  added.  The  derivative  of  the  penalty  function  is  shown  in  (3.16). 
This  equation  was  added  to  the  original  two-norm  gradient  when  the  stability  constraint  was 
violated.  Therefore, 


—  =  2  *  A^  *  max(real(eig(A2)))  *  (— )  (3.16) 

where  ^  is  defined  in  Section  3. 1.2.5. 

3.1 .4  Use  ofFMIN for  Hoo  Norm  Calculation. 

Problem  Definition:  The  highest  percentage  of  computer  time  used  to  calculate  the  H2/ Hoo 
curve  was  spent  doing  the  frequency  search  over  the  singular  value  plot.  Since  the  infinity- 
norm  is  found  by  a  frequency  search  over  the  singular  value  plot  for  the  maximum  magnitude, 
the  finer  the  search  grid  over  this  plot  the  better  the  infinity-norm  calculation.  The  finer  the 
grid,  however,  the  more  expensive  the  calculation  of  the  Hoo  norm. 

Problem  Solution:  Use  of  a  function  in  MATLAB  called  FMIN  to  find  the  maximum  point 
within  defined  bands  for  each  peak  improves  the  run  time.  A  frequency  search  is  still  required 
to  identify  the  peaks  and  set  up  the  bands  for  the  problem,  but  this  search  is  only  done  one 
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time  per  7  step.  The  original  code  did  a  frequency  search  for  every  function  evaluation.  The 
frequency  search  for  the  new  code  can  also  use  many  fewer  points  since  it  only  needs  slope 
information  to  identify  the  peaks  and  generate  the  bands,  but  is  no  longer  required  to  find  the 
precise  value  of  the  infinity-norm. 

Each  function  evaluation  requested  by  the  SQP  optimization  routine  solves  for  the 
Hoo  norm  for  each  band  using  FMIN.  This  function  finds  the  minimum  of  an  unconstrained 
function  of  one  variable  within  a  fixed  interval.  The  function  that  is  to  be  minimized  must 
be  continuous.  Also,  FMIN  may  return  a  local  minimum  if  the  function  it  is  given  is  not 
unimodal.  This  method  was  used  to  find  the  minimum  of  the  negative  of  the  singular  value 
plot  since  FMIN  finds  the  minimum  of  a  function.  It  was  found  that  using  FMIN  provided 
significantly  fewer  evaluations  and  took  less  time  to  run  than  the  previous  frequency  search. 
An  example  of  the  savings  is  where  previously  800  points  were  used  to  search  over  the  singular 
value  plot  to  find  the  maximum  value,  FMIN  would  find  the  maximum  in  approximately  20 
function  evaluations.  This  equates  to  a  97  percent  savings  in  time.  This  comparison  was  made 
when  there  was  only  one  peak  in  the  singular  value  plot,  and  therefore  only  one  band.  FMIN 
also  offers  the  advantage  of  specifying  the  tolerance  to  which  the  infinity-norm  is  calculated 
without  a  significant  increase  in  the  number  of  function  evaluations. 

3.1.5  Seating. 

Problem  Definition:  Numerical  difficulties  are  often  encountered  in  an  optimization  process 
when  one  constraint  is  larger  or  changes  more  quickly  than  another  constraint.  Difficulties  are 
also  encountered  when  the  value  of  the  objective  is  not  of  the  same  order  of  magnitude  as  the 
constraints.  The  goal  then  is  to  normalize  the  objective  and  constraints  to  unity  before  each 
optimization  run  (i.e.,  for  each  7  step  for  the  mixed  H2IH00  problem). 

The  original  code  for  this  calculation  used  scale  factors  to  multiply  the  objective  and 
constraints.  These  scale  factors  were  arbitrarily  set,  depending  on  the  problem  and  the  point 
on  the  H2I  Hoo  curve  where  the  optimization  is  begun. 
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Problem  Solution:  The  objective  was  normalized  by  the  last  value  of  the  two-norm.  A 
constant  scale  factor  was  used  to  multiply  the  normalized  H2  objective.  See  (3.17).  It  was 
found  that  by  scaling  the  objective  constantly  in  this  manner,  the  time  spent  in  the  line  search 
was  much  less.  A  scale  factor  of  (0.001)  was  used  for  all  examples  shown  in  this  thesis.  The 
scale  factor  used  for  the  infinity-norm  was  found  by  normalizing  the  infinity-norm  constraint 
using  the  current  set  constraint  value  of  7.  See  (3.18).  Finally,  a  constant  scale  factor  of 
(1000)  was  used  for  the  stability  constraint.  See  (3.19).  The  stability  constraint  was  not 
normalized.  The  intent  for  this  constraint  was  that  it  always  stay  negative,  not  that  it  approach 
zero.  This  scale  factor,  along  with  the  penalty  function  discussed  earlier,  worked  weU  to 
prevent  a  non-stabilizing  controller  design.  The  scalings  were 


Scalle  for  2-Norm  Objective  = 

(1/1000) 

(3.17) 

Scale  for  Infinity-Norm  Constraint  = 

^/\\Tedg-i  ||oo 

(3.18) 

Scale  for  Stability  Constraint  = 

1000 

(3.19) 

where  ||Tsu,^_j  \\l  and  ||Ted^_j  ||^  are  the  two-norm  calculation  and  the  constraint  level  set  for 
the  infinity-norm  for  the  previous  7  iteration.  The  objective  and  constraints  then  become 


Objective  = 

i/l|r.«,,_J|^(i/iooo)||r,^||^ 

(3.20) 

Infinity-Norm  Constraint  = 

i/\\Tedg_l  ||oo||Ted||oo 

(3.21) 

Stability  Constraint  = 

1000[max  3?e(A)]  <  0 

(3.22) 

3.1 .6  Modifications  for  MIMO  Designs. 

Problem  Definition:  The  code  for  mixed  H2IH00  optimization  needs  to  be  setup  to  calculate 
a  solution  for  a  MIMO  problem.  The  infinity-norm  calculation  must  search  multiple  singular 
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Alt  Singui£U'  Value  Plots 


Figure  3.3  Singular  Value  Plot  for  MIMO  Example  (7=800) 

value  plots  to  find  a  maximum  and  therefore  the  calculation  of  peaks  and  bands  must  be 
adjusted. 


Problem  Solution:  The  modification  that  needed  to  be  made  for  this  case  was  for  the  Hoo  norm 
calculation.  The  maximum  singular  value  over  all  the  singular  value  plots  for  the  MIMO 
system  was  found.  The  FMIN  function  was  still  used  to  calculate  the  peak  values  within 
the  bands.  Two  examples  of  MIMO  singular  value  plots  with  bands  and  peaks  identified  are 
shown  in  Figures  3.3  and  3.4. 

3.1.7  Automatic  Hoo  Steps. 

Problem  Definition:  The  improvements  to  the  code  significantly  reduced  the  run  time  required 
to  generate  the  H2 /Hoo  curve.  However,  large  step  sizes  in  7  caused  numerical  iU  conditioning 
of  the  problem.  This  caused  the  user  to  set  a  step  size,  run  the  optimization  code,  stop  the  code 
when  ill  conditioning  occurred,  save  the  controller  and  restart  the  optimizer  with  a  different 
step  size. 
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All  Singular  Value  Plots 


Figure  3.4  Singular  Value  Plot  for  MIMO  Example  (7=360) 


Problem  Solution:  A  routine  was  added  that  would  automatically  set  the  step  size,  and  therefore 
generate  the  entire  mixed  curve  automatically.  The  autostep  function  begins  by  taking  two  7 
steps  from  H2  optimal.  The  information  from  these  two  steps  is  then  used  to  form  a  linear 
approximation  of  the  mixed  curve.  An  estimate  is  made  using  that  approximation  of  what  size 
the  next  7  step  should  be  for  a  specified  change  in  the  two-norm.  See  Figure  3.5.  If  the  step 
size  is  larger  than  the  original  step  size  it  is  ignored.  If  the  step  size  is  smaller  than  the  original 
7  step  size,  the  next  iteration  is  set  at  that  7  level.  The  original  specified  levels  of  percentage 
changes  of  infinity-norm  and  two-norm  can  be  input  by  the  designer. 

This  automatic  routine  works  well  when  the  designer  chooses  a  small  percentage  change 
in  7  for  the  first  step  size.  The  larger  the  percentage  the  better  the  chance  of  numerical  ill 
conditioning.  The  examples  run  using  this  routine  used  a  change  in  a  of  5  percent  and  a 
change  in  7  of  5  percent. 

3.1 .8  Requiring  H^o  Constraint  to  be  Active. 

Problem  Definition:  A  theorem  was  given  by  Walker  [Wal94]  which  proved  that  if  the  con¬ 
troller  order  is  greater  than  or  equal  to  the  order  of  the  if  2  problem,  then  the  mixed  H2  /iioo 
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Figure  3.5  Linear  Projection  to  Estimate  Gamma  Step 

solution  lies  on  the  boundary  of  the  Hoo  constraint.  It  was  found  that  by  setting  the  opti¬ 
mization  problem  up  with  an  inequality  constraint  occasionally  would  result  in  cases  where 
the  ifoo  constraint  would  not  be  active.  If  the  problem  was  set  up  with  equality  constraints, 
the  optimizer  would  try  to  force  all  the  H(x>  constraints  from  the  frequency  bands  to  zero. 
However,  when  multiple  peaks  were  present,  it  is  not  possible  to  know  in  advance  which 
peak  should  be  constrained  as  an  equality. 

Problem  Solution:  The  mixed  H2I Hoo  active  constraint  requirement  was  that  only  the  maxi¬ 
mum  Hoo  constraint  (from  the  bands)  must  be  active.  Therefore,  a  modification  was  made  to 
the  MATLAB  SQP  routine.  The  modification  required  that  the  maximum  Hoo  constraint  was 
close  to  zero  within  some  tolerance.  This  change  resulted  in  the  Hoo  constraint  being  active  at 
each  step  along  the  mixed  H^  /  Hoo  curve.  It  should  be  noted  that  this  requirement  is  only  true 
for  one  Hoo  constraint.  If  multiple  Hoo  constraints  are  added  to  the  problem,  (i.e.,  UTedi  ||oo 
and  ||red2||oo,  not  constraints  from  several  peaks  of  one  Hoo  problem),  then  only  one  Hoo 
constraint  is  required  to  be  active.  Therefore,  in  the  multiple  Hoo  constraint  case,  this  active 
requirement  is  dropped. 
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3 .2  Numerical  Test  Example 

A  SISO  example  was  chosen  to  demonstrate  these  fixes.  It  does  not  represent  any 
physical  system.  The  example  is  taken  from  [DS79],  where  it  was  used  to  demonstrate  how 
to  recover  robustness  with  observers.  This  study  used  this  example  to  demonstrate  the  fixes  to 
the  optimization  code  because  of  the  limited  number  of  design  variables  required  to  develop 
the  controller.  It  was  also  chosen  because  it  had  a  low  7  at  optimal,  and  would  require 
less  optimization  time  to  generate  the  mixed  H2IH00  curve.  This  was  important  for  ease  of 
making  comparisons  between  the  previous  and  the  improved  numerical  method. 

3.2.1  Problem  Set-Up.  The  equations  for  the  system  can  be  represented  as 


X  =  Ax  +  Bu  + 


■  0  1  ■ 

’O' 

■  35  ■ 

X  -b 

u  -f- 

-3  -4 

1 

-61_ 

y  =  Cx  +  n 
=  [2  1  ]  x  +  n 

where  the  open-loop  plant  transfer  function  for  this  system  is  given  by 


(3.23) 

(3.24) 

(3.25) 

(3.26) 


(3.27) 


5.2.7 .7  H2  Problem.  The  objective  of  the  77^2  design  is  to  find  an  internally 
stabilizing  controller  which  minimizes  the  response  of  the  system  when  disturbances  and 
noises  are  present.  For  this  example,  the  H2  problem  is  a  standard  LQG  problem.  The 
following  matrices  were  chosen  in  constructing  the  regulator  [DS79]: 


2: 

Qc 


Hx 


H^H  =  80 


4x/5[v^ 
35 


^/35 


1  ]  X 
\/35 

1 


(3.28) 

(3.29) 


\ 
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1 


(3.30) 


Rc  — 

These  choices  produce  an  optimal  control  law  of 

u  =  -KcX  =  -[hQ  10]  a;  (3.31) 

which  produces  the  closed-loop  regulator  poles 

A,/  =  -7.0±i2.0  (3.32) 

An  estimator  (i.e.,  Kalman  Filter)  was  designed  to  reconstruct  the  states.  The  following 
matrices  were  selected  for  the  Kalman  Filter  design  [DS79]: 


Qi  = 

Rj  =1 


35 

-61 


(i)[35  -61 ; 


(3.33) 

(3.34) 


The  resulting  H2  matrices  are  as  follows: 


A2 

Bw 


0  1 
-3  -4 
35  O' 

-61  0 


Bu2  — 

a  - 

Cy2  = 


0 

1 

4*a/5 

[2  1] 


1 

0  0 


(3.35) 

(3.36) 

(3.37) 

(3.38) 

(3.39) 
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Dy^  =  [01] 


(3.40) 

(3.41) 


Dyu  —  [0] 


(3.42) 


32.1.2  Hoo  Problem.  The //^oo  problem  is  a  weighted  sensitivity  minimiza¬ 
tion  problem.  It  is  used  to  improve  tracking  performance.  The  closed  loop  sensitivity  was 
chosen  to  give  good  noise  rejection,  tracking  at  low  frequency  and  attenuate  system  response 
to  high  frequency  noise  and  uncertainties.  The  i^oo  problem  can  be  set-up  to  minimize  the 
weighted  sensitivity  of  a  system.  The  sensitivity  weight  was  selected  as  the  inverse  of  the 
desired  sensitivity,  given  by 


W,= 


s  +  12.5 
S-H0.7 


The  state  space  representation  of  the  weight  can  be  written  as 


(3.43) 


A.  = 

Bs  = 
Cs  = 

D,  = 

The  resulting  Hoo  matrices  are  as  follows: 


Bd  = 


1-0.7] 

(3.44) 

11] 

(3.45) 

111.0] 

(3.46) 

[11 

(3.47) 

Bs  *  Cy2 


0 


0 

A., 


(3.48) 

(3.49) 
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Buoo 

= 

Bu2 

0 

(3.50) 

C, 

=  {D,*Cy2  C,] 

(3.51) 

r 

'^yoo 

=  [Cy2 

0] 

(3.52) 

Ded 

=  P.] 

(3.53) 

Dqu 

=  [0] 

(3.54) 

Dyd 

=  [1] 

(3.55) 

=  [0] 

(3.56) 

where  the  submatrices  above  were  given  previously. 

3.2.2  Test  Example  Results.  The  numerical  code,  as  it  was  prior  to  this  work, 
was  run  for  this  example.  The  fuU  mixed  H2/H00  curve  could  not  be  generated.  Due  to 
numerical  problems  discussed  earlier  in  this  chapter,  the  lowest  7  level  that  could  be  achieved 
was  7=4.  The  7  for  this  problem  was  found  using  a  routine  in  MATLAB  called  HINFSYN. 
The  minimum  7  achievable  for  the  third  order  Hoo  problem  was  found  to  be  7=1 .0.  The  mixed 
H2/H00  controller  was  run  at  the  order  of  the  H2  problem,  which  is  2.  Thus,  7=1.0  is  not 
necessarily  achievable  at  order  2,  but  as  will  be  seen,  it  can  be  reduced  to  much  lower  than  4. 

The  problem  of  the  multiple  singular  value  peaks  was  found  to  be  the  reason  the  solution 
close  to  7  could  not  be  reached.  When  the  frequency  search  was  divided  into  two  regions,  (see 
Figure  3.6),  and  separate  -ffco  constraints  were  set  for  each  region,  7=1.001  was  achieved.  A 
comparison  of  the  mixed  curves  before  the  division  for  the  peaks  and  after  is  shown  in  Figure 
3.7  and  Figure  3.8.  The  step  size  in  7  was  unity. 

The  solution  time  using  the  method  of  manually  dividing  the  peaks  and  running  the 
origmal  code  was  approximately  14  hours.  This  is  an  excessive  amount  of  time  when  consid¬ 
ering  that  this  was  a  second  order  control  problem  that  only  required  seven  design  variables 
and  had  a  maximum  of  two  peaks  (i.e.,  two  Hoo  constraints).  Also,  the  number  of  controllers 
generated  in  this  time  frame  was  only  15.  Finally,  this  method  required  a  manual  division 
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Singular  Value  Plot 


Figure  3 .6  Singular  Value  Plot  for  Test  Example 
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Figure  3.7  Mixed  H2/H00  Curve  Original  Code 

of  peaks  by  the  designer.  The  division  of  peaks  was  relatively  simple  for  this  problem  since 
there  were  only  two  peaks,  but  for  problems  that  had  four  or  five  peaks  this  manual  division 
technique  would  be  very  time  consuming. 

Another  comparison  between  the  old  and  new  method  is  the  maximum  number  of 
function  and  gradient  evaluations  required  (see  Table  3.1).  The  old  method,  when  the  peaks 
were  divided  manually,  required  an  average  number  of  function  evaluations  of  approximately 
400  per  7  iteration.  The  new  method  required  only  an  average  of  approximately  10  function 
evaluations  per  7  iteration. 
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Figure  3.8  Mixed  H2IH00  Curve  Improved  Code 


Table  3.1  Number  of  Function  and  Gradient  Evaluations  Comparison 


Original  Code  New  Code 


7 

NUMFUN 

NUMGRAD 

NUMFUN 

NUM  GRAD 

14 

88 

32 

4 

4 

13 

285 

40 

4 

4 

12 

685 

65 

4 

4 

11 

344 

76 

4 

4 

10 

92 

25 

4 

4 

9 

685 

64 

11 

11 

8 

689 

57 

4 

4 

7 

10 

10 

5 

5 

6 

11 

11 

5 

5 

5 

11 

11 

5 

5 

4 

11 

11 

5 

5 

3 

688 

58 

6 

6 

2 

12 

12 

6 

6 

1 

694 

231 

119 

98 

TOTAL 

4305 

703 

186 

165 
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Figure  3.9  Mixed  H2IH00  curve  with  over  150  controllers 


Note  that  the  maximum  number  of  function  evaluations  per  7  step  was  set  at  700.  There 

were  five  cases  in  the  original  code  where  number  of  function  evaluations  was  close  to  700; 

/ 

in  these  cases,  the  optimizer  was  stopped  at  maximum  iterations.  The  value  of  the  size  of 
the  line  search  step  for  each  of  these  cases  was  consistently  negative,  indicating  the  solution 
to  within  the  specified  tolerance  would  not  have  been  found.  Controllers  with  infinity-norms 
very  close  to  tolerance  were  generated  in  this  manner. 

Using  the  improvements  discussed  in  Section  3.1,  the  run  time  was  reduced  from  14 
hours  down  to  2  minutes  as  measured  in  CPU  time.  The  CPU  time  compared  very  closely  to 
real  time  when  working  on  a  dedicated  SPARC20  SUN  workstation.  Many  more  controllers 
could  also  be  calculated.  The  time  required  to  calculate  the  150  controllers  shown  in  Figure 
3.9  was  17  CPU  minutes. 

33  Summary 

This  chapter  presented  eight  general  improvements  made  to  the  optimization  code  for 
the  mixed  H2III00  controller  design.  It  discussed  how  each  problem  was  recognized  and  a 
solution  implemented.  These  improvements  were  then  tested  on  a  small  example  to  compare 
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improved  efficiency.  The  example  shows  that  a  solution  for  the  mixed  problem  could  be 
found  very  close  to  7,  where  previously  this  controller  could  not  be  found.  The  mixed 
H2  / Hoo  controller  came  very  close  to  the  absolute  minmum  7,  but  could  not  equal  it  since  it 
was  computed  at  an  order  less  than  the  order  of  the  Hoo  problem.  The  example  also  showed 
that  the  new  code  ran  approximately  200  times  faster  than  the  previous  code.  This  was 
calculated  for  the  example  where  the  frequency  search  was  manually  divided  and  a  complete 
mixed  curve  could  be  generated.  Finally,  this  example  demonstrates  the  speed  in  which  many 
controllers  along  the  mixed  H2/ Hoo  curve  could  be  found.  The  results  showed  it  took  17 
CPU  minutes  for  150  controllers  to  be  calculated  with  the  new  code  and  14  CPU  hours  for  15 
controllers  to  be  calculated  with  the  old  code  and  the  manually  divided  peaks. 
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IV.  F-16  Example 


4.1  Introduction 

The  purpose  of  this  chapter  is  to  design  mixed  H2IH00  controllers  for  a  normal  accel¬ 
eration  command  following  model  of  the  AFTI  F-16.  This  example  was  first  examined  by 
Luke  [Luk93].  He  designed  both  LQG/LTR  and  mixed  H2/H00  controllers  for  this  problem. 
Luke’s  study  was  done  to  compare  the  performance  of  the  two  controllers  and  draw  conclu¬ 
sions  about  the  effectiveness  of  the  mixed  H2IH00  design.  After  completion,  it  was  found 
that  the  bandwidth  for  the  LTR  design  was  different  than  the  bandwidth  for  the  sensitivity 
weighting  of  the  Hoo  design.  The  maximum  bandwidth  for  the  LQG/LTR  design  was  10 
rad/sec,  whereas  the  Hoo  sensitivity  weighting’s  bandwidth  was  1  rad/sec.  The  objective  for 
this  thesis  was  to  change  the  Hoo  sensitivity  weighting’s  bandwidth  to  10  rad/sec,  rerun  the 
mixed  H2I Hoc  controllers,  and  compare  them  to  the  LQG/LTR  design  calculated  by  Luke. 

Another  objective  for  this  example  is  to  examine  controllers  of  orders  different  than 
that  of  the  H2  problem.  Three  different  order  controllers  were  evaluated.  The  first  was  at  the 
order  of  the  H2  problem  and  the  second  was  the  order  of  the  Hoo  problem.  Finally,  a  reduced 
order  mixed  H2I  Hoo  controller,  which  was  one  order  lower  than  the  H2  problem,  was  also 
evaluated. 

This  chapter  first  defines  the  F-16  problem  set-up  in  Section  4.2.  The  mixed  H2I Hoo 
controllers  for  different  orders  are  defined  and  analyzed  in  Section  4.3.  Finally,  the  compar¬ 
isons  between  the  LQG/LTR  and  the  mixed  H2I Hoo  designs  are  presented  in  Section  4.4. 

4.2  Problem  Set-Up 

This  example  is  a  normal  acceleration  command  following  model  of  an  AFTI  F-16 
flying  at  sea  level,  Mach  0.6.  The  linear  model  consists  of  a  servo,  longitudinal  dynamics 
and  a  delay.  This  is  shown  in  Figure  4.1.  A  more  detailed  model  of  the  F-16  model  is  shown 
in  Figure  4.2  This  section  will  first  discuss  the  state  space  development  for  the  H2  and  Hoo 
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Figure  4,2  F-16Plant 


problems.  Next,  the  noise  model  used  for  this  example  will  be  developed.  Finally,  a  brief 
description  of  an  LQG/LTR  design  will  be  presented. 

The  states  for  the  core  plant  are 

u 

a 

(4.1) 

0 

q 


In  (4.1),  the  states,  u,  a,  0  and,  q  are  perturbations  in  the  forward  speed,  angle  of  attack,  pitch 
angle,  and  pitch  rate,  respectively.  All  angles  and  angular  rates  are  in  radians  and  radians! sec 
and  velocity  is  in  ft! sec. 
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4.2.1  Servo.  The  servo  is  modeled  as  a  first  order  lag  with  a  pole  at  -20.  This  is 


represented  in  transfer  function  form  as 


^  _  20 
^ec  S  -|-  20 

This  can  be  written  in  state  space  format  as 


(4.2) 


xs,  =  [-20]xi.  -I-  [20] 


(4.3) 


Se  =  [ijXi,  [0]6e, 


(4.4) 


4.2.2  Longitudinal  Dynamics.  The  longitudinal  dynamics  of  the  F-16  at  sea  level 
flying  at  Mach  0.6  were  obtained  from  Baird  [Bai92]. 


-0.0015 

37.382 

-32.200 

-17.940 

0.0021 

-0.0001 

-1.4910 

-0.0013 

0.9960 

^cp 

-0.1880 

0 

0 

0 

1.0000 

0 

-0.0004 

9.7530 

0.0003 

-0.9600 

-19.0400 

(4.5) 


0.0015  35.2640  0.0272  -0.3340 


Xcp  +  [— 4.3660]5e 


(4.6) 


4.2  J  Time  Delay.  For  the  time  delay,  a  period  of  0.05  seconds  is  chosen,  which 
leads  to  a  first  order  Padd  approximation  transfer  function  of 


N,  _  40 -g 
■^zg  5  -f-  40 


(4.7) 


The  corresponding  state  space  representation  of  the  delay  is 


id  =  [-40]a:d  -|-  [IjiV^^ 

(4.8) 

N,  =  [80]xd  4-  [-IjAT,^ 

(4.9) 
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4.2.4  Complete  State  Space.  The  servo,  longitudinal  dynamics,  and  delay  can  be 
combined  into  one  complete  state  space.  The  states  Xcp  represent  those  for  the  core  plant 
(the  longitudinal  dynamics).  The  states  xs^  represents  those  for  the  servo,  and  the  states  Xi 
represent  those  for  the  time  delay. 


Xcp 

A 

■^cp 

Bcp  Cs^ 

0 

Xcp 

BcpDs^ 

XSe 

= 

0 

ASe 

0 

Xs, 

+ 

Bs. 

Xd 

Bd^cp 

Bd  Bcp  B 

Xd 

BdDcpDs^ 

(4.10) 


DdCcp  DdDcpCs^  Cd 


Xcp 

Xs, 

Xd 


+  [DdDcpDs^]6, 


Sc 


(4.11) 


or,  substituting  in  the  values. 


Xp  — 


N, 


'  -0.0015 

37.3820 
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0.0021 

0 

’  0  ' 

-0.0001 

-1.4910 

-0.0013 

0.9960 

-0.1880 

0 

0 

0 

0 

0 

1.0000 

0 

0 

Xp~\~ 

0 

-0.0004 

9.7530 

0.0003 

-0.9600 

-19.040 

0 

0 

0 

0 

0 

0 

-20 

0 

20 

0.0015 

35.2640 

0.0272 

-0.3340 

-4.366 

-40.000 

0 

^del 


-0.0015  -35.264  -0.0272  0.3340  4.3360  80.000 


(4.12) 

Xp  +  [0]8cc  (4.13) 


4 .2 .5  Noise.  The  noise  models  for  the  LQG/LTR  design  are  given  in  this  subsection. 
These  noises  were  used  in  all  simulations.  The  noises  introduced  in  this  system  were  wind 
gusts  and  sensor  noise.  The  wind  gusts  were  introduced  as  perturbations  in  angle  of  attack. 
They  were  modeled  as  UWGN  inputs  that  were  then  filtered.  The  sensor  noise,  which  entered 
at  the  output  of  the  plant,  was  also  modeled  as  a  filtered  UWGN. 
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The  filter  used  for  the  wind  gust  was  developed  by  Baird  [Bai92].  It  is  a  low  pass  filter 


which  cuts  off  at  6.7  rad/sec.  The  filter  in  transfer  function  form  is 


^  .  w  0.0187 

^wsisoi^)  ~  ~  ^  g  y  (4-14) 

The  sensor  noise  filter  was  designed  to  closely  match  the  assumed  error  in  measuring 
a  1  G  response.  This  assumed  error  was  ±  0.04Gs  or  4  percent.  The  filter  is  high  pass  with 
scaling  to  adjust  the  rms  to  0.04.  The  filter  in  transfer  function  form  is 


n  0.04(s  +  0.01) 

~  10{s  +  10) 


(4.15) 


4.2.6  LQG/LTR  Design.  A  brief  discussion  of  Luke’s  [Luk93]  LQG/LTR  design 
is  given  here.  This  background  is  presented  so  that  the  design  process  for  LQG/LTR  is 
understood  when  comparing  the  mixed  H2IH00  controller  to  the  LQG/LTR  controller.  For 
details  on  this  design  see  [Luk93]. 

The  first  step  in  the  design  process  is  to  develop  a  desired  loop  shape.  A  loop  shape 
was  chosen  that  had  a  maximum  bandwidth  of  10  radians/sec  and  a  20  dB/decade  roll-off  at 
bandwidth.  In  order  to  make  the  loop  shape  meet  the  specified  barriers,  a  bank  of  integrators 
is  added  to  the  design  plant.  These  integrators  increase  the  order  of  the  controller  by  two. 

The  next  step  was  to  design  the  Kalman  filter  and  the  regulator.  The  design  parameters 
H,  p,  A  and  p  (which  are  the  state  weighting,  the  control  weighting,  the  process  noise  distri¬ 
bution  matrix,  and  the  sensor  noise  intensity,  respectively)  were  chosen  by  Luke  [Luk93].The 
design  variable  H  was  set  equal  to  the  C  matrix  of  the  plant,  />=1,  A  was  set  equal  to  the 
second  column  in  the  A  matrix,  and  p  was  chosen  to  be  100.  The  variable  that  affects  recovery, 
q,  was  varied,  with  q  =  1000  chosen.  The  resulting  controller  will  be  compared  with  mixed 
H2IH00  controller  designs  in  Section  4.4. 
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4 .3  Mixed  H2  / Hoo  Controller  Design 

The  mixed  H2lH^  controller  design  method  allows  the  designer  to  choose  levels  of 
tracking  and  noise  rejection  performance.  The  objective  for  this  example  is  to  calculate  mixed 
H2IH00  curves  for  different  orders  of  the  controller,  and  then  compare  their  performance. 
The  mixed  controller  designs  at  different  orders  will  also  be  compared  to  the  LQG/LTR 
design  from  the  previous  section.  The  mixed  design  is  developed  using  the  updated  numerical 
optimization  code  discussed  in  Chapter  3.  This  section  is  broken  into  five  parts.  First,  the  H2 
and  Hoo  problems  are  developed.  Then  different  controller  orders  are  analyzed  for  the  mixed 
H2IH00  design.  The  first  controller  order  evaluated  is  the  order  of  the  H2  plant.  Next  the 
mixed  H2  /  Hoo  controller  with  order  of  the  Hoo  plant  is  examined,  and  finally,  a  reduced  order 
controller  one  order  less  than  the  H2  problem  is  analyzed. 


4.3.1  H2  Problem.  The  minimum  order  of  the  mixed  H2  / Hoo  controller  for  which 
strong  statements  can  be  made  about  its  characteristics  is  driven  by  the  order  of  the  H2  plant. 
Therefore,  the  designer  may  wish  to  keep  the  order  of  the  H2  problem  as  low  as  possible.  To 
keep  the  order  as  low  as  possible,  it  is  best  to  use  static  weights  on  z  and  w.  When  static 
weights  are  used,  the  H2  problem  then  becomes  the  standard  LQG  problem.  This  means  the 
filters  used  for  the  noise  will  attenuate  all  frequencies  equally.  The  values  chosen  for  the  static 
gains  were  (5o=0.0028  for  the  process  noise  (i.e.  wind)  and  i?/=0.004  for  the  sensor  noise. 
The  control  weighting  chosen  was  i2c=100  and  H  was  chosen  equal  to  the  C  matrix  of  the 
plant.  The  vector  F  corresponds  to  the  column  in  the  plant  A  matrix  that  multiplies  the  a 
state.  The  value  A  is  equal  to  1. 

The  H2  problem  for  this  example  can  then  be  written  as 
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The  full  state  space  matrices  are  in  Appendix  A. 

43.2  Hoo  Problem .  The  tracking  performance  of  the  controller  for  a  step  input 

is  based  on  the  weighting  put  on  the  system  sensitivity.  Usually  this  weighting  is  chosen  to  be 
the  inverse  of  the  desired  sensitivity.  Since  a  bandwidth  of  close  to  10  rad/sec  was  required  to 
compare  to  the  LQG/TTR  design,  the  weighting  chosen  was 


W,{s) 


s  +  8.0001 
s  +  0.0001 


(4.16) 


Figure  4.3  shows  that  the  loop  shape  meets  the  boundaries  chosen  for  the  target  loop  shape. 
In  state  space  the  weight  can  be  written  as 


A,  =  [-0.0001]  Ba  =  [1]  C,  =  [8]  Ds  =  [1] 


Now,  with  the  sensitivity  weight  added  to  the  Hoo  system,  the  resulting  complete  Hoo 
state  space  matrices  can  be  written  as 
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The  full  state  space  matrices  for  this  example  are  in  Appendix  A.  Since  the  sensitivity  weight 
added  to  this  problem  is  not  strictly  proper,  the  Ded  term  is  not  zero.  Also,  Deu=0,  which 
means  the  problem  is  singular.  However,  from  Walker’s  development  [Wal94]  there  is  no 
requirement  for  Ded  to  be  zero  and  singular  Hoo  problems  are  allowed. 

433  Hi  Order  Mixed  Hi / H^  Solution.  The  optimization  routine  was  initiated  at 
the  Hi  optimal  design.  The  calculation  of  the  Hi  order  mixed  controller  previously  required 
days  to  compute  only  six  points  on  the  mixed  Hi  (Hoo  curve.  The  new  code  generated  over 
150  controllers  in  less  than  3  hours.  A  plot  of  the  mixed  curve  is  shown  in  Figure  4.4. 

Note  that  the  lowest  value  of  7  achieved  for  this  example  was  1.361.  The  minimum 
value  of  7  found  for  the  Hoo  problem,  which  is  one  order  higher  than  the  Hi  order,  was 
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Sensitivity  F-16  PhugoW 


Figure  4.5  Sensitivity  of  Closed  Loop  System  :  Dashed  Line  -  Target  Loop  Shape  (Order 
H2) 

7=1.3201.  This  value  should  not  be  expected  to  be  achieved,  since  this  is  for  a  controller  of 
the  Hcc  order.  As  the  order  of  the  controller  is  decreased,  7  increases. 

Sensitivity  curves  for  this  problem  are  plotted  in  Figure  4.5.  This  shows  the  trade-off 
between  H2  and  /foo-  As  7  is  decreased,  the  if 00  part  of  the  problem  becomes  more  dominant 
and  the  sensitivity  plot  approaches  the  target  loop  shape,  (i.e.,  the  dashed  line  on  Figure  4.5). 
The  margins  for  the  problem  can  be  found  using  this  plot.  Notice  the  peak  is  relatively  low 
indicating  good  margins. 

4.3.4  Hoo  Order  Mixed  H2  / ifoo  Solution.  Since  the  optimal  value  of  7  at  the  ffoo 
order  of  the  plant  is  known  to  be  1.3201,  the  optimizer  at  this  order  should  get  very  close  to 
this  value.  It  was  desired  to  look  at  the  effect  of  performance  of  a  mixed  controller  as  close 
as  possible  to  this  optimal  7  level.  The  optimization  code  must  be  started  with  a  controller  of 
the  appropriate  order.  The  previously  generated  controllers  have  order  one  less  than  required 
here.  The  order  of  the  controller  can  be  increased  to  provide  an  acceptable  starting  point  by 
adding  a  near  pole-zero  cancellation  or  by  adding  a  pole  beyond  the  bandwidth  of  the  system. 
The  second  approach  was  chosen.  A  lag  was  added  to  a  controller  found  for  the  H2  order 
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Singular  Value  Plot 


Figure  4.6  Singular  Value  Plot  During  Optimization  (Order  Hoo) 

problem  with  a  7=1.5.  Consequently,  the  lag  increased  7  for  the  starting  controller  from  7=1.5 
to  7=1.58.  This  lag  was 


1000 
s  + 1000 


(4.17) 


Running  the  optimization  code  at  the  order  of  the  Hoo  problem  should  allow  the  optimizer  to 
get  very  close  to  the  optimal  7  level.  The  final  7  level  achieved  for  this  problem  was  1.331. 
Figure  4.6  shows  the  resulting  singular  value  plot  of  the  closed  loop  system  for  7=1.331. 
There  were  six  peaks  constrained  for  that  singular  value  plot.  It  should  be  noted  how  flat  the 
singular  value  plot  has  become.  A  plot  of  the  mixed  H2I Hoo  curve  at  the  Hoo  order  for  the 
controller  is  shown  in  Figure  4.7. 


The  sensitivity  plot  for  this  order  controller  shows  that  the  low  frequency  dip  is  removed 
as  7  is  decreased.  It  also  shows  the  higher  frequency  dip  in  the  sensitivity  is  not  removed. 
This  plot  is  given  in  Figure  4.8. 


43.5  Reduced  Order  H^ /Hoo  Solution.  The  final  question  is  what  is  the  controller 
performance  when  the  order  of  the  controller  is  reduced  from  H2  optimal.  A  direct  technique 
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Figure  4,7  Mixed  Curve  (Order  Hoo) 


Sensitivity  F-16  Phugoid 


Figure  4.8  Sensitivity  of  Closed  Loop  System  (Order  Hoo) 

was  used  to  reduce  the  order  of  the  controller.  This  was  accomplished  by  running  the 
MATLAB  function  ’SCHMR’.  This  function  uses  the  Schur  method  for  model  reduction. 
The  function  SCHMR  constrains  the  infinity-norm  of  the  difference  between  the  original  and 
reduced  transfer  functions  to  be  less  than  a  specified  value.  The  designer  also  has  the  option  of 
choosing  the  order  of  the  reduced  transfer  function.  There  is  no  guarantee  that  the  controller  is 


4-12 


1.^ 


1.8h  Initial  Controller 


1.1  i - ' - ' - ' - ' - * 

0  500  1000  1500  2000  2500 

Gamma 

Figure  4.9  Mixed  H2  / Hoo  Ciu^'e  (Order  Reduced) 

stabilizing  after  this  model  reduction  routine  has  been  performed.  The  reduced  order  controller 
was  formed  from  a  mixed  H2  /  Hoo  design  that  was  the  order  of  the  H2  problem.  This  controller 
had  a  7=2(K)0  and  a=1.156.  This  controller  was  found  to  be  both  stable  and  stabilizing  for 
the  closed  loop  system.  The  mixed  H2I  Hoo  curve  for  the  reduced  order  problem  is  shown  in 
Figure  4.9.  The  point  on  the  graph  marked  with  an  ’0’  was  the  starting  point  for  the  reduced 
order  controller.  At  this  point,  the  two-norm  was  1.725  and  the  infinity-norm  was  2000.  Note 
the  infinity-norm  did  not  change  from  the  sixth  order  controller  to  the  fifth  order.  The  optimum 
two-norm  value  was  found  for  this  controller  by  setting  the  7  step  level  equal  to  zero  and 
letting  the  optimizer  find  the  controller  with  the  lowest  two-norm.  The  optimum  two-norm 
for  the  starting  controller  is  marked  by  a  on  Figure  4.9  at  a  value  of  1.156. 

The  H2  optimal  controller  at  orders  less  than  the  plant  is  no  longer  unique.  There  may 
be  many  different  controllers  at  this  order  that  give  the  same  H2  norm,  but  may  have  different 
Hoo  norms  for  T^d-  No  attempt  was  made  to  find  an  optimal  reduced  order  H2  controller  here; 
rather,  the  starting  point  described  above  was  used. 

A  plot  of  the  closed  loop  sensitivity  is  shown  in  Figure  4.10.  The  sensitivity  plot  for 
reduced  order  shows  a  slight  improvement  from  the  H2  order  in  being  able  to  match  the  target 
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Figure  4. 10  Sensitivity  of  Closed  Loop  System  (Order  Reduced) 

loop  shape.  The  phugoid  dip  in  the  sensitivity  plot  is  removed  completely  at  the  lowest  7 
controller  for  the  reduced  design.  The  order  sensitivity  plot  shows  the  short  period  dip 
stUl  remains  at  the  lowest  7  level  achieved. 

4  A  Comparison  of  Controller  Designs 

Each  of  the  mixed  H2IH00  designs  was  compared  to  the  LQG/LTR  design  discussed 
in  Section  4.2.6.  Simulations  were  run  on  MATLAB  SIMULINK.  Comparisons  were  made 
for  each  controller’s  closed-loop  response  to  a  unit  step  and  an  initial  5  degree  perturbation  in 
angle  of  attack. 

4A.1  Mixed  H2  Order  vs  LQGILTR .  This  section  presents  the  comparison  of  the 
LQGA-TR  controller  versus  the  H2  order  mixed  H2IH00  controller.  Figure  4.11  shows  the 
mixed  curve  for  the  H2  order  controller.  The  LQG  controller  is  shown  with  the  symbol  ’0’ 
in  Figure  4.11.  Note  that  the  LQG  two-norm  is  higher  than  that  of  H2‘,  this  is  due  to  the 
augmented  integrators  in  the  LQG  design.  Since  Loop  Transfer  Recover  techniques  were 
used  to  improve  the  tracking  of  the  LQG  design,  the  LQG/LTR  design  is  shown  on  the  mixed 
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Figure  4. 1 1  Mix  H2 / Hoo  Curve  with  LQG  design  (Order  H2) 


Gamma 


Figure  4.12  Mix  H2IH00  Curve  with  LQGA-TR  design  (Order  H2) 

H2IH00  curve  in  Figure  4.12.  It  is  important  to  note  that  the  LQG  and  LQG/LTR  designs  are 
eighth  order  where  the  mixed  design  shown  here  is  sixth  order. 

Figure  4.12  shows  that  the  eighth  order  LQG/LTR  controller  can  be  improved  in  both 
noise  performance  and  tracking  compared  to  the  mixed  6th  order  controller.  The  LQG/LTR 
design  resulted  in  a=1.58  and  7=551.88.  A  mixed  controller  was  chosen  which  resulted  in 
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Table  4.1  Gain  and  Phase  Margin  Comparison  (H2  Order) 


Mixed  H2I  Hoo 

VGMs  (dB) 

VPMs  (deg) 

VGMt(dB) 

VPMt(deg) 

Total  VGM(dB) 

Total  VPM  (deg) 

[-4.4  9.8] 

±39.5 

[-10.3  4.5] 

±40.8 

[-10.3  9.8] 

±40.8 

LQG/LTR 

VGMs  (dB) 

VPMs  (deg) 

VGMt  (dB) 

VPMt(deg) 

Total  VGM(dB) 

Total  VPM  (deg) 

[-4.6  10.6] 

±41.35 

[-6.5  3.7] 

±30.8 

[-6.5  10.6] 

±41.35 

a=1.76  and  7=1.5  to  compare  to  the  LQG/LTR  design.  The  gain  and  phase  margins  were 
calculated  for  the  mixed  controller  and  the  LQG/LTR  controller,  and  are  shown  in 

Table  4.1.  The  total  phase  margins  are  virtually  the  same  for  both  controllers,  but  the  mixed 
H2IH00  total  gain  margins  are  slightly  better  than  LQG/LTR. 

Figure  4.13  shows  the  response  of  both  controllers  for  an  initial  5  degree  perturbation 
in  angle  of  attack  and  a  1  G  unit  step  in  N^.  The  noise  characteristics  look  about  the  same  but 
the  mixed  H2IH00  controller  is  a  better  tracker.  The  mixed  controller  takes  less  time  to  settle 
for  the  response  to  an  initial  five  degree  angle  of  attack.  It  also  shows  less  of  an  overshoot  for 
the  step-response. 


4.4.2  Mixed  Hoo  Order  vs  LQGILTR.  The  mixed  curve  at  the  order  of  the 

Hoo  problem  is  shown  in  Figure  4.14.  The  point  ’0’  on  the  curve  is  the  LQG/LTR  8th  order 
design.  As  in  the  previous  example,  there  is  a  mixed  design  at  7th  order  that  improves  the 
two-norm  and  infinity-norm  of  the  system  compared  to  the  LQG/LTR  controller.  The  mixed 
controller  chosen  for  comparison  produced  a  two-norm  equal  to  1.66  and  an  infinity-norm 
equal  to  1.5.  The  gain  and  phase  margins  were  calculated  for  the  mixed  H2I  Hoo  controller. 
These  are  in  Table  4.2.  The  margins  for  Hoo  order  are  approximately  the  same  as  the  margins 
for  H2  order. 

The  response  of  the  mixed  design  at  H^o  order  is  shown  in  Figure  4.15.  The  response 
of  the  mixed  design  demonstrates  better  performance  than  the  LQG/LTR  design.  The  mixed 
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Figure  4.13  Responses  of  Mixed  H2 1 Hoo  (Order  H2)  (dash)  vs  LQG/LTR  (solid) 


Table  4.2  Gain  and  Phase  Margin  Comparison  Order) 


Mixed  H2I  Hoo 

VGMs(dB) 

VPMs(deg) 

VGMt(dB) 

VPMt(deg) 

Total  VGM(dB) 

Total  VPM  (deg) 

[-4.4  9.6] 

±39.2 

[-10.3  4.5] 

±40.7 

[-10.3  9.6] 

±40.7 

LQG/LTR 

VGMs  (dB) 

VPMs  (deg) 

VGMt  (dB) 

VPMt  (deg) 

Total  VGM(dB) 

Total  VPM  (deg) 

[-4.6  10.6] 

±  41.35 

[-6.5  3.7] 

±30.8 

[-6.5  10.6] 

41.35 
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Figure  4. 14  Mix  H2  / Hao  Curve  with  LQG/LTR  design  (Order  Hoo) 

controller,  as  in  the  H2  order,  has  a  quicker  settling  time  for  the  initial  condition  and  less  of  an 
overshoot  for  the  step  response.  The  noise  performance  is  approximately  the  same.  Overall, 
there  is  no  real  benefit  in  using  the  H^o  order  versus  the  H2  order  mixed  controller. 

4.43  Mixed  Reduced  Order  vs  LQGILTR  .  The  mixed  00  curve  for  one  order 
lower  than  the  H2  problem  (order  5)  is  shown  in  Figure  4.16.  The  mixed  controller  chosen  for 
comparison  gave  a  two-norm  of  1.89  and  an  infinity-norm  of  1.59.  The  LQG/LTR  design  is 
also  plotted  on  this  curve  for  comparison  purposes.  It  can  be  seen  that  even  at  this  lower  order 
for  the  mixed  controller  the  two-norm  and  infinity-norm  can  be  improved  over  the  LQG/LTR 
design.  The  gain  and  phase  margins  were  also  calculated  for  this  mixed  reduced  order  H2 /Hoo 
controller.  They  are  in  Table  4.3.  Both  the  gain  and  phase  margins  are  slightly  better  for  the 
LQG/LTR  design  than  the  mixed  H2I  Hoo  design.  The  lower  order  mixed  margins  are  slightly 
worse  than  the  H2  and  i/00  order  margins  with  little  degredation  in  tracking  performance. 
This  was  expected  at  reduced  order. 

The  response  performance  for  the  reduced  order  controller  versus  the  LQG/LTR  design 
is  shown  in  Figure  4.17.  This  shows  that  the  mixed  design  has  a  faster  response  time  to  an 
initial  input  of  a  5  degree  alpha  command  and  the  unit  step  response  to  a  1  G  input  of  Nz  has 
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Figure  4.15  Responses  of  Mixed  H2I Hoo  (Order  Hoo)  (dash)  vs  LQG/LTR  (solid) 
Table  4.3  Gain  and  Phase  Margin  Comparison  (5th  Order) 


Mixed  H2IH00 

VGMs  (dB) 

VPMs  (deg) 

VGMt(dB) 

VPMt  (deg) 

Total  VGM(dB) 

Total  VPM  (deg) 

[-4.3  9.23] 

±38.2 

[-9.9  4.5] 

±39.8 

[-9.9  4.5] 

±39.8 

LQG/LTR 

VGMs  (dB) 

VPMs  (deg) 

VGMt(dB) 

VPMt  (deg) 

Total  VGM(dB) 

Total  VPM  (deg) 

[-4.6  10.6] 

±  41.35 

[-6.5  3.7] 

±30.8 

[-6.5  10.6] 

±  41.35 
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Figure  4. 1 6  Mix  H2  /  Hoc  Curve  with  LQG  design  (Reduced  Order) 

less  overshoot.  This  mixed  H2IH00  design  of  order  five  performs  better  than  the  LQG/LTR 
design  of  order  eight. 


4.4.4  Summary  of  Comparisons.  The  purpose  of  this  example  was  to  determine 
the  effectiveness  of  mixed  H2lH^  optimization  versus  LQG/LTR  for  a  real-world  example. 
The  results  from  this  example  show  the  mixed  approach  is  better  than  LQG/LTR  design 
method.  This  example  demonstrated  that  the  mixed  H2IH00  controller  at  three  orders  less 
than  the  LQG/LTR  design  is  a  better  tracker.  The  noise  response  for  the  different  order  mixed 
controllers  versus  the  LQG/LTR  remained  relatively  the  same.  The  optimal  infinity-norm  at 
seventh  order  was  1.33,  sixth  order  was  1.36  and  fifth  order  was  1.59.  This  similar  noise  and 
tracking  response  between  the  different  order  controllers  is  due  to  choosing  mixed  controllers 
with  the  same  approximate  two-norm  and  infinity-norm.  Although  this  was  not  done  by 
design,  the  infinity-norm  for  the  mixed  controllers  was  1.5, 1.5  and  1.59.  The  two-norm  for 
the  designs  chosen  were  1.66, 1.76  and  1.89  for  seventh,  sixth  and  fifth  order  mixed  H2IH00 
controllers.  Each  choice  was  based  on  achieving  the  best  two-norm  and  infinity-norm  for  that 
particular  curve. 
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Figure  4.17  Responses  of  Mixed  H2IH00  (Reduced  Order)  (dash)  vs  LQG/LTR  (solid) 


The  results  presented  in  this  chapter  also  demonstrate  the  improvement  in  the  H^fHoo 
optimization  code.  This  optimization  improvement  provides  the  capability  to  trade  off  H2  and 
Hoo  objectives  for  the  entire  mixed  H2IH00  curve.  Previously  the  mixed  solution  could  not 
generate  enough  of  the  curve  in  the  region  of  the  optimal  mixed  solution  to  examine  the 
full  tradeoff. 

The  results  of  this  example  showed  improvement  in  performance  at  lower  order  con¬ 
trollers.  The  results  also  demonstrated  the  improved  numeric  capability  to  generate  the  fuU 
mixed  H2IH00  curve  in  significantly  less  time  than  was  previously  required.  These  were 
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important  for  a  SISO  example;  however,  these  improvements  become  much  more  significant 
for  a  MIMO  example.  The  missile  MIMO  problem  is  addressed  in  Chapter  V. 
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V.  MIMO  Example 


The  control  problem  discussed  in  this  chapter  is  for  an  air-to-air  missile  that  uses  tail- 
control  fins.  The  purpose  of  this  chapter  is  to  design  a  mixed  H2IH00  controller  for  this 
example  and  then  compare  that  design  to  a  LQG/LTR  design  previously  completed  by  Brown 
[Bro91].  The  mixed  H2/ Hoc  design  was  attempted  with  the  previous  optimization  code,  but 
due  to  numerical  difficulties  no  results  could  be  generated.  The  improved  optimization  code 
was  applied  to  this  MIMO  problem  to  test  its  effectiveness.  The  full  H2IH00  curves  could  be 
generated  for  this  problem.  Since  the  goal  is  to  make  a  comparison  between  the  LQG/LTR 
and  the  mixed  controller  designs,  this  example  concentrates  on  one  point  in  the  missile  flight 
envelope. 

5.1  Background 

The  added  challenge  of  the  controller  design  for  this  example  is  that  the  plant  is  non¬ 
minimum  phase  with  lightly  damped  complex  poles.  Although  this  is  usually  a  difficult 
problem  for  LQG/LTR,  Brown  successfully  calculated  the  LQG/LTR  controller  using  the 
method  of  eigenstructure  reassignment  through  static  feedback.  This  static  feedback  matrix 
does  not  add  to  the  dimension  of  the  plant,  therefore  keeping  the  controller  order  at  a  minimum. 
Eigenstructure  reassignment  was  used  first  on  the  plant  to  shift  the  poles  and  then  LQG/LTR 
was  applied  to  the  problem.  This  method  proved  beneficial  to  the  design. 

5.1.1  Eigenstructure  reassignment.  It  has  been  shown  by  Andry  [ASC83]  that 
complete  pole  reassignment  is  possible  if  the  system  is  completely  observable,  completely 
controllable,  andn<r-fm-l  where  n  is  the  order  of  the  system,  r  is  the  number  of  inputs, 
and  m  is  the  number  of  outputs.  This  missile  MIMO  example  did  not  require  all  the  poles 
to  be  moved.  A  method  developed  in  [Bro91]  allows  a  limited  number  of  the  system’s  poles 
to  be  reassigned  with  the  remaining  poles  drifting  to  unassigned  positions.  This  method  also 
allows  arbitrary  assignment  of  portions  of  the  system’s  eigenvectors. 
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The  method  begins  with  the  state  space  representation  of  a  completely  controllable 


completely  observable  system. 


X  ==  -4x  +  5u 
y  =  Cx 


(5.1) 


where  A  is  of  rank  n,  B  is  rank  m,  and  C  is  rank  r.  Note  this  development  assumes  that  there 
is  no  D  matrix. 


The  problem  then  is  to  calculate  a  constant  gain  matrix  Kstat  that  has  dimensions  mxr, 
which  will  reassign  the  poles  through  inner  loop  feedback.  The  state  space  equation  with 
Kstat  is  now 


X  —  X  +  ;Bu 

y  =  Cx 


(5.2) 


5.1.2  Eigenvector  assignability.  The  definition  of  eigenvector  assignability  was 
addressed  by  Andry  [ASC83].  He  discusses  the  extent  of  possible  eigenvector  specification 
for  the  closed-loop  system.  The  definition  of  the  eigenvalue,  eigenvector  pair  (A,  and  v, )  is 
given  as: 

[A  +  BKstatC]vi  =  X{y/i  (5.3) 

or,  rearranged, 

v.  =  [XJ  -  A]~^  BKstatCvi  (5.4) 


It  is  required,  from  the  development,  that  the  eivenvector  v,-  must  lie  in  the  subspace 
spanned  by  the  columns  of  [A,/  —  A]~^B  which  is  rank  B.  Since  the  desired  eigenvector  will 
probably  not  be  in  the  exact  subspace  defined  an  achievable  eigenvector  is  calculated  by  a 
projection  of  the  desired  vector  onto  the  subspace  spanned  by  the  columns  of  [Xil  —  A]B. 

Only  some  of  the  components  of  the  eigenvector  are  reassignable.  The  desired 
assignable  components  are  chosen  by  the  designer  to  try  to  affect  performance  goals.  The  de- 
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sired  eigenvector  can  be  reordered  to  contain  the  assignable  components  at  the  top  of  the  vector 
and  the  unassignable  components  at  the  bottom  of  the  vector.  This  allows  the  calculation  of 
the  assignable  eigenvector. 

5.13  Kstat  Calculation.  The  gain  matrix  Katat  is  developed  assuming  the  eigen¬ 
vectors  are  in  the  assignable  general  subspace.  The  process  of  forming  Katat  begins  by 
transforming  the  system  so  that  the  B  matrix  becomes: 

[Im] 

...  (5.5) 

.  . 

A  similarity  transformation  matrix,  T,  is  formed  such  that 

X  =  Tx'  (5.6) 

where 

T  =  [  [S]  i  [P]  ]  (5.7) 

This  similarity  transformation  allows  the  eigenvalues  of  the  system  to  be  unchanged  but  the 
eigenvectors  have  been  transformed  by: 

v'  =  T-\  (5.8) 

The  equation  for  the  closed-loop  eigenvalues  and  eigenvectors  written  in  partitioned  form 
becomes 

[Xilm  ~  *411]  [”*^12]  [^i]  [.^m]  [Vi] 

—  ^stat^  p.y) 

[-3^21]  [Xiln-m  -  A22]  J  [  K]  J  [  [0]  J  [  [Wi]  _ 
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Using  the  top  equation,  rearranging  and  defining  Ai 
written  as 


All  '•  Ai2 


the  equation  can  be 


[Ai  Vj  —  ^if]i 


(5.10) 


This  equation  must  be  true  for  all  eigenvalue/assignable  eigenvector  pairs.  Rearranging  (5.10) 
and  solving  for  Kstat  gives  the  equation 


=  [N  -  AW]  [CW]-' 


(5.11) 


5.2  Missile  Model 

The  missile  is  assumed  rigid  and  the  missile’s  airframe  dynamics  are  modeled  using  six 
state  variables;  u,  a,  p,  q,  and  r.  These  variables  represent  forward  velocity  in  feet/second, 
angle  of  attack  in  the  missile’s  pitch-plane  in  radians,  sideslip  angle  in  radians,  roll  rate,  pitch 
rate,  and  yaw  rate  in  radians/second.  The  model  is  driven  by  four  fin  inputs:  (5i,  S2, 63,  and  ^4. 
The  outputs  are  p,  q,  r,  N^,  Ny,  and  N^.  The  outputs  p,qmdr  are  measurements  from  the 
rate  gyros  of  roll,  pitch  and  yaw  respectively.  The  terms  N^,  Ny  and  Nz  are  accelerometer 
measurements  from  the  rate  gyros  in  the  x,  y  and  z  directions  repsectively. 

Brown  developed  the  linearized  equations  of  motion  for  this  problem.  The  linearized 
equations  were  compared  to  the  sponsor’s  six  degree  of  freedom  (6-DOF)  model.  Overall, 
it  was  determined  that  the  linearized  model  was  adequate  for  the  LQG/LTR  study.  These 
equations  of  motion  were  put  into  standard  state  space  format  using  the  notation 


X  =  Ax  -f  Bu 

y  =  Cx  -1-  Du  (5.12) 
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Figure  5.1  Schematic  of  Closed-Loop  System  for  Missile 


where 
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The  A,  B  and  C  matrices  for  this  problem  are  shown  in  Appendix  B. 

53  LQGILTR  Design 

The  design  process  includes  evaluating  performance  boundaries,  evaluating  the  need  for 
augmenting  the  plant  with  integrators,  eigenstructure  reassignment  and  the  filter  and  regulator 
design.  The  closed-loop  system  can  be  represented  by  Figure  5.1.  The  inputs  to  Kstat  are  p, 
q,  r,  Ny  and  N^.  The  output  of  Kstat  is  then  fed  directly  into  the  plant  before  the  bank  of 
integrators.  The  performance  bounds  for  this  design  include  a  68  percent  rise  time  in  normal 
acceleration  within  200ms.  Another  requirement  is  minimum  overshoot.  A  final  requirement 
is  zero  steady  state  error  for  a  step  input. 
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The  system  response  equation  was  then  defined  by  Brown.  The  particular  quantities 
that  were  to  be  controlled  were  roll  rate(p),  y-axis  acceleration  (Ny)  and  z-axis  accleration 
(N^).  This  is  a  subset  of  the  system  outputs.  Defining  the  response  equation  in  this  manner 
effectively  discards  the  output  information  of  the  two  states.  The  system  response  equation  is: 


z  =  Hx  = 


P 

Ny 

N, 


(5.13) 


Order  reduction  of  the  problem  was  first  accomplished  by  recognizing  that  the  first  state, 
u  was  totally  decoupled  from  the  other  states  and  did  not  appear  in  any  of  the  equations  for 
the  system  response  inputs.  Therefore,  the  linear  model  could  be  reduced  by  one  state.  This 
was  done  by  eliminating  the  rows  in  the  state  space  A,  B,  and  C  matrices  that  corresponded 
with  the  state  u.  Finally  an  order  reduction  in  the  plant  was  achieved  by  reducing  the  order  of 
the  model  of  the  actuator  dynamics.  The  actuator  dynamics  of  the  missile  can  be  modeled  as 


)fm  ~  ^^Xm  T"  'BynUm 


where 


(5.14) 


(5.15) 


Brown  reduced  the  order  of  the  actuator  models.  The  original  model  presented  by  the 
sponsor  was  second  order  with  high  damping  (C  =  0.8)  and  undamped  natural  frequency 
of  408  rad/sec.  He  showed  that  a  useful  approximation  of  this  actuator  was  a  first  order 
model  with  the  pole  at  s=-400.  The  state  space  representation  of  the  plant  augmented  with  the 
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actuators  is  given  as: 


Xjn 

>Am 

BmCa 

+ 

0 

Xa 

0 

Aa 

_ 

Ym 


V  C 


(5.16) 


The  four  fin  deflections  (^i  through  ^4)  were  expressed  in  standard  fin  deflections 
8g,  Sr).  These  were  the  required  fin  deflections  to  produce  roll-rate,  pitch-rate,  and  yaw- 
rate  respectively.  The  transformation  from  four  inputs  to  three  was  done  by  mixing  the  fin 
deflections  accordingly 


(5.17) 


JJ.J  Design  for  Condition  7.  Flight  condition  7  represents  the  missile  flying  at 
Mach  2  near  sea  level.  Brown  did  this  design  both  with  and  without  eigenstructure  reassign¬ 
ment.  This  discussion  will  only  present  the  design  with  the  eigenstructure  reassignment  since 
that  design  performed  better. 

The  eigenvector  reassignment  was  done  with  a  goal  of  achieving  specific  relationships 
between  the  states.  Reassigning  the  defined  components  for  the  eigenvectors  is  the  same  for 
each  flight  condition.  The  one  real  pole  was  moved  to  the  left  slightly  to  keep  gains  low.  The 
two  complex  pairs  of  poles  were  moved  slightly  to  the  left,  increasing  the  natural  frequency 
and  particularly  the  damping.  Original  poles  (-1.6,  —1.22  ±  i  13. 9  and  -1.31  ±  j  14.0).  The 
new  poles  were  (-2,  — 15  ±  il5.0  and  —15  ±  jl5.0). 
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The  filter  was  designed  with  a  value  of  //=3000  to  accomplish  a  30  dB  lowering  of  the 
singular  value  plot.  The  controller  design  used  g^=10000  to  get  the  best  loop  shape.  This 
design  only  slightly  violated  the  left  boundary  for  the  loop  shape. 

5.4  Mixed  H2IH00  Design 

The  mixed  ^^2/^00  design  was  approached  in  two  ways.  The  first  formulation  follows 
Brown’s  LQG/LTR  design  method  where  two  of  the  five  outputs  are  discarded.  The  two 
measurements  not  used  were  q  and  r.  This  will  be  referred  to  as  the  3x3  MIMO  problem.  The 
second  approach  is  not  to  ignore  the  other  two  measurements  and  include  them  in  the  mixed 
H2IH00  design.  This  will  be  called  the  5x3  MIMO  problem. 

The  LQG/LTR  design  by  Brown  is  14th  order,  where  the  plant  to  be  controlled  is  eighth 
order.  Since  the  mixed  H2  /  design  is  the  same  order  of  the  plant  of  the  H2  problem,  the 
resulting  controller  will  be  order  eight. 

5.4.1  3x3  MIMO  Example.  This  development  ignores  two  measurements.  This 
problem  is  broken  into  the  separate  H2  and  Hoo  problems. 

5.4.1 .1  H2  Problem.  The  goal  of  the  H2  problem  is  to  keep  the  order  of  the 
problem  as  low  as  possible  since  that  will  be  the  order  of  the  controller  for  the  mixed  design. 
The  weights  that  will  be  used  are  static  and  again  as  in  the  SISO  example  the  set-up  becomes 
the  standard  LQG  problem.  The  order  for  the  H2  problem  is  then  eight. 

The  weight  chosen  for  the  LQG  design  was  H  =  Cpianu  so  that  Qc  =  H^H.  The 
control  weighting  Rc  was  equally  weighted  for  each  fin.  This  was  done  due  to  the  fact  that 
the  missile  is  symmetric  about  its  longitudinal  axis.  Thus,  Rc  =  pi  where  /)=5000  for  this 
design. 

The  Kalman  filter  design  used  MIMO  noise  models.  These  models  were  developed  in 
the  same  way  as  the  SISO  noise  model  development  in  Chapter  4.  The  difference  between  the 
wind  noise  model  for  the  SISO  problem  and  the  noise  model  for  the  MIMO  problem  was  that 
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a  perturbation  in  13  must  also  be  added.  The  process  noise  was  modeled  as 


G 


^MIMO 


2x10-3  0 

0  2x10-'' 


(5.18) 


and  no  filter  dynamics  were  added.  The  sensor  noise  was  developed  in  exactly  the  same  way 
as  the  SISO  problem.  The  sensor  noise  model  used  was 
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The  Kalman  filter  design  then  used 

(5.19) 
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The  H2  problem  is  formulated  as 
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(5.20) 

(5.21) 


(5.22) 


(5.23) 
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(5.24) 
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Appendix  B  contains  the  full  state  space  matrices. 


5.4.12  Hoo  Problem.  The  sensitivity  weighting  used  for  the  ifoo  problem 
was  derived  by  the  specifications  used  by  Brown  in  designing  the  LQG/LTR  controller.  This 
ensures  a  better  one-for-one  comparison  between  the  two  methods.  The  target  loop  transfer 
function  chosen  was 

20 

Gis)K{s)  = - (5.25) 

the  sensitivity  weighting  matrix  can  be  written  as 


W,{s)  = 
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(5.26) 


the  state  space  representation  of  this  weighting  matrix  is  given  by 
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The  Hoc  part  of  the  3x3  MIMO  problem  is  then 
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(5.27) 


(5.28) 


(5.29) 


The  order  of  the  Hoo  problem  is  11.  The  values  in  the  state  space  matrices  are  given  in 
Appendix  B. 


5.4.2  5x3  MIMO  Example.  This  set-up  uses  all  five  measurements  of  the  system.  It 
is  very  similar  to  the  3x3  MIMO  set-up.  The  H2  design  is  exactly  the  same  as  the  3x3  MIMO 
problem.  The  same  values  for  H,  Qc,  Rc,  Qo  and  Rf  were  used.  The  only  difference  is  that 
the  C  matrix  for  the  plant  now  has  five  measurements. 

5.4.2.1  Hoo  Problem.  The  most  significant  difference  between  the  5x3 
MIMO  design  and  the  3x3  MIMO  design  is  the  weighting  matrix  used.  The  five  outputs  make 
the  weighting  matrix  5x5,  but  the  requirements  for  tracking  haven’t  changed.  There  is  no 
requirement  to  track  pitch  and  yaw  rates  {q  and  r);  therefore,  the  sensitivity  weights  for  q  and 
r  are  set  to  one.  These  measurements  will  then  be  used  for  maigin  control  only  and  not  to 
improve  the  tracking.  The  weighting  matrix  Wg  is  then 
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this  can  be  represented  in  state  space  format  as 
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The  set-up  for  the  Hoo  problem  is  symbolically  the  same  as  the  3x3  MIMO  problem.  The 


weight  changed  but  the  order  of  the  Hoo  problem  remained  1 1  since  the  two  weights  added 


were  static.  The  state  space  matrices  with  numerical  values  for  this  problem  are  located  in 
Appendix  C. 


5.5  Mixed  H2I Hoo  Controller  Results 

The  optimization  of  this  problem  is  more  difficult  than  the  SISO  problem  due  to  the 
number  of  design  variables.  The  SISO  problem,  at  the  order  of  the  H2  problem,  has  23  design 
variables.  The  MIMO  3x3  case  has  63  design  variables.  The  MIMO  5x3  case  has  79  design 
variables.  The  increase  in  the  number  of  design  variables  from  the  3x3  to  the  5x3  is  caused 
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Figure  5.2  Mixed  H2IH00  Curve  (3x3  MIMO) 

by  the  inaeased  size  of  the  controller  C  matrix  due  to  the  two  additional  measurements.  No 
mixed  H2I if 00  controllers  were  able  to  be  found  using  the  original  optimization  code.  Both  of 
these  problems  have  been  run  using  the  improved  optimization  code  and  full  if2/ifco  curves 
were  generated  for  both  examples.  The  mixed  curves  along  with  the  sensitivity  plots  are 
presented  in  this  section. 

5.5.1  3x3  Mixed  H2  /  ifoo  Controller  Design.  The  mixed  H2  / if<x>  controUer  design 
for  this  problem  was  initialized  at  the  H2  controller.  The  order  of  this  controller  is  eight. 
This  3x3  MIMO  example  is  what  will  be  used  to  compare  to  the  LQG/LTR  design  since 
two  measurements  were  not  used  for  either  design  case.  The  mixed  H2lH^  curve  for  this 
problem  is  shown  in  Figure  5.2.  The  minimum  7  achieved  for  this  problem  was  202.85.  The 
minimum  absolute  7  for  the  problem  with  an  11th  order  controller  is  1.783.  The  minimum  7 
at  8th  order  was  expected  to  be  much  closer  to  1.783.  This  could  be  due  to  two  factors;  1)  the 
optimization  routine  is  stuck  in  a  local  minimum  of  the  design  space  and  can  not  move  away 
from  that  point  or  2)  the  MATLAB  ’CONSTR’  can  not  handle  the  numerics  of  the  problem  at 
that  point  in  the  mixed  H2 1  Hoo  design  space. 
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Figure  5.3  Singular  Value  Plot  (3x3  MIMO) 


The  singular  value  plot  of  the  closed-loop  sensitivity  approaches  the  target  sensitivity 
as  7  approaches  7  for  the  system.  This  is  shown  in  Figure  5.3.  The  controller  with  7=202.85 
produced  approximately  the  same  shape  for  the  sensitivity  plot  as  the  inverse  of  the  weight, 
but  the  plot  was  shifted  up  by  about  30  dB.  This  indicates  the  tracking  performance  will  be  less 
than  desired.  Also,  higher  values  of  7  produced  sensitivity  plots  closer  to  the  target  sensitivity. 

5.5.2  5x3  Mixed  H2 Controller  Design.  The  mixed  H2  /  H^o  controller  design 
for  this  problem  was  started  at  the  H2  controller.  The  order  of  this  controller  is  eight.  This 
design  contains  all  five  outputs.  The  mixed  H2 /Hoo  curve  for  this  problem  is  shown  in  Figure 
5.4.  The  minimum  7  achieved  for  this  problem  was  346. 

The  maximum  singular  value  of  the  sensitivity  approaches  the  target  sensitivity  as  7 
nears  the  optimal  value  for  the  system.  This  is  shown  in  Figure  5.5.  The  controller  with  an 
7  equal  to  346  produced  the  sensitivity  curve  that  has  the  same  shape  as  the  target  sensitivity, 
but  is  shifted  up  by  about  30  dB. 
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Figure  5.4  Mixed //2/jf^oo  Curve  (5x3  MIMO) 
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Figure  5.5  Singular  Value  Plot  (5x3  MIMO) 


5.6  Comparison  of  Mixed  H2  /  controller  designs  to  LQGILTR 

The  mixed  H2  / Hoo  controller  with  a  7=202.85  was  simulated  for  a  lateral  acceleration 
command  of  1  G.  The  simulation  showed  the  tracking  performance  to  be  very  poor.  A 
simulation  of  the  three  outputs  for  a  1  G  step  in  Ny  is  shown  in  Figure  5.6.  The  LQG/LTR 


5-15 


Time  (sec) 


Figure  5.6  Response  of  system  to  1  G  iVj,  step  (Mixed  H^jHoo  Controller) 
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Figure  5.7  Response  of  system  to  1  G  iVj,  step  (LQG/LTR  Controller) 

controller  with  Kstat  design  for  the  same  step  input  produced  very  good  tracking,  as  shown  in 
Figure  5.7. 

5.6.1  Comparison  of  Mixed  H2  /  Hoo  with  Kstat  to  LQGILTR.  The  mixed  H2  /  Hoo 

design  was  done  for  the  original  missile  plant  with  the  non-minimum  phase  zeros  and  the 
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Figure  5.8  Response  of  system  with  Kstat  to  a  1  G  step  (Hoo  Optimal  Controller) 

lightly  damped  poles.  This  is  a  very  hard  plant  to  control.  Brown  used  the  static  gain  matrix 
to  move  the  poles  of  the  system.  Brown’s  study  showed  that  without  moving  the  poles,  the 
LQG/LTR  controller  did  not  perform  well.  Therefore,  the  next  step  was  to  check  the 
optimal  controller  design  for  the  modified  plant,  (i.e.,  with  the  static  gain  closed  in  a  loop  with 
the  plant).  The  simulation  of  the  controller  showed  very  good  tracking  performance.  This  is 
shown  in  Figure  5.8.  The  optimal  controller  order,  however,  was  eleven.  The  H2  optimal 
order  for  the  problem  was  eight. 

The  mixed  H2I  Hoo  design  was  then  done  for  the  modified  missile  plant  using  Kstau 
and  still  for  order  eight.  The  mixed  H2I Hoo  curve  is  shown  in  Figure  5.9.  It  was  simulated 
for  the  same  1  G  step  acceleration  in  Ny.  The  tracking  performance  was  about  the  same  as 
the  design  for  the  original  missile  plant.  There  is  an  unacceptable  amount  of  overshoot  and 
extremely  long  settling  time  for  the  response  to  a  1  G  step.  A  plot  of  the  response  is  in 
Figure  5.10. 

5.6.2  Comparison  of  Mixed  H2  / Hoo  at  11th  order  to  LQGILTR.  Previous  results 
showed  that  the  Hoo  optimal  controller  performed  well.  See  Figure  5.8.  The  final  attempt 
made  to  design  an  H2I  Hoo  controller  was  to  run  the  optimization  code  at  1 1th  order  (i.e.,  the 
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Figure  5.9  Mixed  H2IH00  Curve  with  Kstat  (3x3  MIMO) 
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Figure  5.10  Response  of  Mixed  H^jHoo  Controller  to  IG  Ny  step  with  Kstat  (H2/H0C, 
Controller) 
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order  of  the  Hoo  problem).  Since  using  Katat  to  do  the  mixed  H2IH00  design  at  8th  order 
did  not  improve  tracking  performance,  Katat  was  not  incorporated  into  the  1 1th  order  design. 
At  11th  order,  7=1.76.  Therefore,  it  is  known  that  the  optimization  routine  should  be  able  to 
get  very  close  to  7=1.76  at  this  order.  The  H2  optimal  mixed  controller  for  the  3x3  problem 
was  increased  from  8th  to  1 1th  order  by  adding  a  lag  to  each  transfer  function.  The  lag  can  be 
represented  as 

-1000  0  0 

Alag  =  0  -1000  0 

0  0  -1000 

1000  0  0 

Clag  =  0  1000  0 

0  0  1000 

0  0  0 
Dlag  =  0  0  0 

0  0  0 

The  mixed  H2/H00  curve  at  11th  order  requires  87  design  variables.  This  higher  order 
problem  controller  design  was  attempted  using  the  MATLAB  ’CONSTR’  optimization  rou¬ 
tine.  Although  the  previous  3x3  and  5x3  problems  at  H2  order  (8th  order)  could  be  run  using 
the  MATLAB  optimizer,  the  3x3  problem  at  11th  order  could  not  be  done  using  MATLAB. 
No  mixed  H2IH00  controllers  could  be  found  using  MATLAB.  Next,  a  FORTRAN/EMSL 
optimization  routine  was  used  for  this  problem.  A  FORTRAN  program  was  written  that 
called  the  same  MATLAB  files  used  for  the  previous  examples  only  where  the  MATLAB 
program  had  called  the  ’CONSTR’  optimization  routine  the  FORTRAN  program  now  called 
IMSL/DNCONG.  For  further  explanation  on  the  details  of  MATLAB  versus  IMSL’s  opti¬ 
mization  routines,  see  Chapter  8.  This  routine  was  able  to  find  mixed  H2IH00  designs  for 
this  problem.  The  mixed  H2fHoo  curve  is  shown  in  Figure  5.11.  The  local  minimum  in 
the  plot  could  not  be  removed  by  the  means  discussed  in  Chapter  6.  The  lowest  7  found 
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Figure  5.12  Singular  Value  Plot  Mixed  H2  / H^o  1 1th  Order  (7=66 

using  the  optimization  routine  was  7=66.  The  absolute  7  for  this  problem  was  1.76.  The 
difference  between  these  two  numbers  was  much  less  than  for  the  8th  order  mixed  design.  The 
closed-loop  singular  value  plot  at  7=66  is  shown  in  Figure  5.12.  Note  how  flat  the  singular 
value  plot  has  become. 
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Figure  5.13  Comparison  of  LQG/LTR  Controller  to  Mixed  H2  / Hoo  Controller 

A  comparison  was  then  made  between  the  a  and  7  for  the  LQG/LTR  design  by  Brown 
with  Katat  to  that  found  for  the  mixed  H2IH00  controller  with  Kstat  at  11th  order  with  the 
lowest  7  found.  This  is  shown  in  Figure  5.13.  This  comparison  shows  two  things:  1)  that  the 
mixed  H2/H00  design  can  never  achieve  the  a  and  7  produced  by  the  LQG/LTR  controller 
or  2)  that  the  mixed  H2IH00  controller  at  7=66  is  at  a  local  minimum  and  could  be  reduced 
further.  The  latter  is  believed  to  be  the  case. 

The  step  response  for  the  mixed  H2IH00  controller  to  a  1  G  step  in  Ny  improved  as 
shown  in  Figure  5.15.  However,  when  compared  with  the  LQG/LTR  design,  the  mixed  has 
five  times  as  much  overshoot  as  LQG/LTR. 

5.7  Summary 

The  mixed  H2IH00  controller  design  was  completed  for  both  the  3x3  MIMO  and  the 
5x3  MIMO  8th  order  problems  using  the  MATLAB  optimizer  ’CONSTR’.  The  optimization 
code  was  able  to  find  mixed  controllers  for  this  larger  MIMO  example.  Mixed  controllers  can 
be  found  in  approximately  a  day  for  this  size  of  problem(79  design  variables).  The  1 1th  order 
could  not  be  optimized  using  MATLAB  ’CONSTR’.  Rather,  a  FORTRAN/IMSL  optimizer 


B«st  Mixed  11th  Order  Controller 
without  Kstat 


LQG/LTR  Controller  with  Kstat 
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Figure  5.14  Response  of  system  to  1  G  TVj,  step  (Mixed  H2  jHoo  Controller  Order  11) 
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Figures. 15  Response  of  system  to  1  GiVj,step(Mixedif2/^^oo  Controller  Order  11)  versus 
LQG/LTR 
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was  used  to  generate  the  mixed  H2  / Hoo  curve.  Due  to  the  size  of  the  problem,  the  1 1th  order 
problem  took  several  days  to  converge  to  7=66. 

Simulation  of  these  controllers  shows  that  the  controller  does  not  track  a  step  input  very 
well  at  8th  order.  There  is  an  unacceptable  amount  of  overshoot  and  the  settling  time  does  not 
meet  the  requirements.  When  the  optimal  order  Hoo  controller  was  tested  on  this  problem, 
it  performed  very  well,  but  it  was  three  orders  higher  than  the  H2  problem.  The  11th  order 
mixed  H2/ Hoo  controller  at  7=66  demonstrated  improved  tracking  performance  compared 
to  8th  order  but  was  still  not  as  good  as  LQG/LTR.  This  is  believed  to  be  a  local  minimum 
problem.  Work  is  currently  being  done  to  help  the  numerical  code  "kick  out  of  local  minima. 
Finally,  it  was  shown  that  moving  the  poles  of  the  system  did  not  improve  the  mixed  H2  /Hoo 
controller  performance  as  it  did  for  the  LQG/LTR  design. 
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VI.  Evidence  of  Local  Minimum 


6.1  Introduction 

It  is  usually  very  difficult,  when  applying  optimization  techniques,  to  be  sure  a  global 
minimum  has  been  found.  This  is  true  for  any  problem  that  is  not  convex.  The  necessary 
conditions  and  a  positive  definite  Hessian  matrix  ensure  that  the  design  found  by  the  opti¬ 
mization  routine  is  a  local  minimum.  If  the  problem  is  not  convex,  the  only  way  to  ensure 
that  the  design  is  a  global  optimum  is  to  prove  that  the  Hessian  is  positive  definite  for  all 
possible  values  of  design  variables.  This  is  rarely  practical  in  a  design  application.  Posing  the 
optimization  problem  as  a  convex  problem  (if  possible)  is  the  best  way  to  ensure  a  design  is  a 
global  minimum. 

The  optimal-order  mixed  H2/H00  controller  design  was  posed  as  a  convex  problem 
by  Walker  [Wal94].  He  states  that  if  the  transfer  function  Tz^  is  in  the  space  H2  then  the 
two-norm  of  Tzw{,Q)  is  a  strictly  convex  functional  of  Q.  If  a  function  is  convex  there  is  one 
global  minimum  and  therefore  the  design  at  that  point  is  unique.  The  optimal  order  for  the 
mixed  problem  has  not  be  determined  analytically.  Numerical  results  have  suggested  that 
the  optimal  order  is  larger  than  the  combined  H2  and  plant.  WeDs  and  Ridgely  [WR92] 
have  suggested  that  the  optimal  order  is  infinite.  The  optimal-order  mixed  H2  /  H^o  controller 
design  is  a  convex  problem  in  the  design  space  of  the  Youla  parameter  Q,  but  what  about  the 
fixed-order  controller  problem?  The  two-norm  objective  and  the  infinity-norm  and  stability 
constraints  are  not  convex  functions  in  the  design  variable  space.  Currently,  it  is  not  know 
whether  it  is  possible  to  pose  the  fixed-order  mixed  problem  in  a  convex  setting. 

This  section  present  results  for  both  the  F-16  SISO  problem  and  the  missile  MIMO 
problem  that  show  evidence  that  the  fixed  order  problem  is  not  convex  as  posed  in  this  work. 
There  are  mixed  H2IH00  curves  generated  for  both  these  examples  which  show  that  local 
minimum  occur. 
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Figure  6. 1  Mixed  /H^o  Curve  Local  Minimum  {H2  Order) 

6.2  Examples  of  Local  Minimum 

Both  the  F-16  and  missile  examples  will  be  discussed  in  this  section.  The  local  minimum 
are  more  frequent  and  more  pronounced  in  the  missile  example. 

6.2.1  F-16  Example.  TheF-16example  was  discussed  in  Chapter  IV.  This  example 
was  the  first  to  show  evidence  of  a  local  minimum.  The  local  minima  for  this  example  occur 
more  frequently  when  small  steps  in  7  are  taken.  They  also  occur  more  frequently  the  closer 
the  design  gets  to  H^o  optimal.  Figure  6.1  shows  a  place  in  the  mixed  H2IH00  curve  where  a 
local  minimum  occurred. 

Further  evidence  that  this  point  was  a  local  minimum  was  the  singular  value  plot.  At 
the  point  the  local  minimum  occurred,  the  shape  of  the  singular  value  plot  changed  drastically. 
This  lead  to  the  conclusion  that  a  completely  different  controller  design  had  been  found,  and 
that  a  local  minimum  was  the  reason  for  this  design  change.  Figure  6.2  shows  the  change  in 
the  singular  value  plot. 

6.2.2  Missile  Example.  The  missile  MEMO  example  also  demonstrated  local 
minima.  The  trends  are  similar  to  the  F-16  case,  in  that  they  occur  more  frequently  with  small 
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Singular  Value  Plot 


Figure  6.2  Singular  Value  Plot  where  local  minimum  occurred  Order) 


Figure  6.3  Mixed  H2  / H^o  Curve  Local  Minimum  (3x3) 


steps  in  7.  They  also  occur  more  frequently  the  closer  the  mixed  design  gets  to  Hoo  optimal. 
A  plot  with  local  minima  for  the  MEMO  3x3  example  is  shown  in  Figure  6.3. 

The  MEMO  5x3  example  also  showed  local  minimum.  A  plot  of  the  mixed  H2IH00 
curve  for  this  example  is  in  Figure  6.4.  An  example  of  local  minimum  occurring  very  close 
to  optimal  is  shown  in  Figure  6.5. 
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Figure  6.4  Mixed  H2 /Hoo  Curve  Local  Minimum  (3x5) 


Figure  6.5  Mixed  H2 1  Hoo  Curve  Local  Minimum  (3x5) 


6.3  Solutions  to  Local  Minimum 

One  way  to  prevent  local  minima  from  occurring  is  by  taking  large  steps  in  7,  thereby 
possibly  skipping  over  a  local  minimum.  This  is  not  a  good  solution,  however,  since  large 
steps  may  cause  numerical  ill  conditioning.  Smaller  steps  work  well  to  make  sure  low  points 
in  the  design  space  are  found.  The  drawback  of  small  steps  is  that  there  is  more  of  a  chance  the 
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optimizer  will  be  stuck  in  a  local  minimum  and  will  not  be  able  to  move  enough  to  find  a  lower 
local  minimum  or  the  global  minimum.  Methods  such  as  simulated  annealing  offer  the  ability 
to  move  away  from  a  local  minimum  and  possibly  find  a  global  minimum.  Although  they 
were  not  used  for  this  study,  they  could  possibly  decrease  the  chances  of  the  design  getting 
stuck  in  a  local  minimum  point  when  taking  small  steps. 

Once  a  local  minimum  has  been  found,  there  is  a  way  to  "fix"  the  mixed  H2 / if 00  design 
curve.  This  is  done  by  using  the  controller  that  was  found  just  after  the  drop  in  the  mixed 
H2I Hoc  curve.  Using  this  controller,  the  drop  in  the  mixed  H2IH00  curve  can  be  "removed" 
by  tracing  the  curve  in  the  opposite  direction,  (i.e.,  taking  positive  steps  in  7),  This  method 
was  used  to  remove  the  local  minimum  points  for  the  F-16  and  missile  examples.  This  is  the 
reason  the  examples  in  Chapters  4  and  5  do  not  show  the  local  minimum  points  as  described 
in  this  chapter. 

6.4  .  Summary 

Local  minimum  have  occurred  in  both  the  F-16  and  missile  example  completed  for 
this  thesis.  The  reason  these  local  minimum  have  not  been  obvious  before  is  due  to  the 
limited  reliability  of  the  code.  When  only  seven  to  ten  controllers  could  be  calculated  for 
the  entire  mixed  H2IH00  curve,  the  chances  of  seeing  this  problem  were  very  low.  With  the 
improved  reliability  and  the  means  to  generate  hundreds  of  controllers  for  a  specific  example, 
the  problem  has  become  more  evident. 

These  examples  demonstrated  that  the  chances  of  a  local  minimum  becoming  obvious 
are  increased  when  smaller  step  sizes  in  7  are  taken.  Local  minima  also  appear  more  frequently 
the  closer  the  design  gets  to  if 00  optimal.  There  is  no  means  to  determine  if  the  curve  generated 
is  the  absolute  lowest  local  minimum.  However,  if  a  local  minimum  is  obvious  to  the  designer, 
(i.e.,  the  mixed  curve  drops  to  a  lower  part  of  the  curve),  that  local  minimum  point  can  be 
removed  from  the  mixed  curve.  This  clearly  shows  that  the  fixed  order  problem  is  not  convex 
in  the  design  variable  space. 
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Vll.  Time  Savings  Comparisons 

The  main  objective  of  this  thesis  was  to  reduce  the  time  to  compute  a  mixed  H2lH^ 
controller.  To  verify  that  this  was  done,  a  time  comparison  was  made  between  the  old  software 
and  the  new  improved  software  for  two  examples.  This  was  done  for  the  small  SISO  system 
discussed  in  Chapter  3  and  for  the  F-16  problem.  No  time  comparison  could  be  made  for  the 
missile  problem,  since  no  solution  could  be  reached  with  the  original  code. 

7.1  SISO  Second  Order  System 

Table  3.1  showed  the  number  of  required  function  evaluations.  This  section  wiU  present 
a  comparison  of  the  required  amount  of  CPU  time  in  seconds  for  the  example  to  be  run  on  a 
dedicated  SPARC20  SUN  workstation  using  the  old  and  new  optimization  code. 

Table  7.1  shows  the  second  order  SISO  example  discussed  in  Chapter  3  with  the  old 
code  and  the  peaks  manually  divided.  This  comparison  was  chosen  since  a  complete  curve 
could  not  be  generated  without  the  peaks  constrained  separately.  Note,  as  discussed  in  Chapter 
3,  the  original  code  with  separately  constrained  peaks  had  five  out  of  the  fourteen  cases  where 
the  iteration  terminated  by  maximum  iterations  exceeded.  This  means  the  tolerances  set  on 
the  design  variables,  the  function  and  the  constraints  were  not  met.  Therefore,  the  code  would 
have  taken  even  longer  to  run  if  the  limit  on  maximum  number  of  iterations  was  not  imposed. 

7.2  F-16  Example 

The  original  code  was  run  for  the  F-16  example.  Few  points  were  generated  due  to  the 
time  required  to  generate  one  point.  A  comparison  of  the  time  to  generate  those  points  using 
the  old  code  versus  the  new  code  is  shown  in  Table  7.2.  To  get  an  answer,  the  scale  factors  of 
the  original  code  were  set  at  Ai=0.001,  A2=l  and  A3=1000.  For  other  scale  factors,  such  as 
0.01, 1  and  10  respectively,  a  controller  for  the  first  7  step  could  not  be  found. 
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Table  7.1 


Time  Comparison  of  Old  versus  New  Optimization  Code  (Test  Example) 


OLD  CODE 

NEW  CODE 

7 

CPU  sec 

CPU  sec 

14 

436.36 

8.16 

13 

1795.41 

3.39 

12 

2130.45 

3.38 

11 

4895.68 

3.58 

10 

4835.21 

3.69 

9 

4903.26 

6.8 

8 

1949.76 

3.41 

7 

1974.93 

4.25 

6 

167.21 

3.98 

5 

521.26 

4.07 

4 

4816.35 

4.05 

3 

464.28 

4.71 

2 

5224.44 

4.72 

1 

3238.2 

85.58 

TOTAL  37352.8  143.71 


Table  7.2  Time  Comparison  of  Old  versus  New  Optimization  Code  (F-16  Example) 


7 

OLD  CODE 
CPU  sec 

NEW  CODE 
CPU  sec 

76744 

518.10 

17.21 

75744 

31.45 

11.91 

74755 

91.80 

14.01 

TOTAL 

641.35 

43.13 

7.5  Summary 

It  is  obvious  that  improving  the  gradients  alone  for  this  problem  would  increase  the 
reliability  and  efficiency  of  the  code.  The  tables  in  this  chapter  give  a  perspective  of  just 
how  much  time  is  saved  from  the  original  method  compared  to  the  new  optimization  code. 
Examples  which  originally  took  hours  now  take  minutes. 
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VIll.  MATLAB  Compared  to  FORTRAN  IMSLIDNCONG 


The  purpose  of  this  chapter  is  to  compare  MATLAB  to  FORTRAN  IMSL/DNCONG 
for  mixed  H2lH^  optimization.  Both  MATLAB  and  the  FORTRAN  IMSL  routines  use 
SQP  as  the  method  of  optimization.  However,  the  SQP  method  is  implemented  differently 
in  each  routine.  First,  a  background  on  IMSL  and  its  implementation  of  SQP  is  given. 
Next,  a  background  on  the  MATLAB  optimization  package  and  its  implementation  of  SQP 
is  presented.  Following  the  background  of  these  methods,  a  comparison  between  the  SQP 
implementations  is  made.  Finally,  the  F-16  example  is  used  to  compare  run  times  and 
convergence  between  MATLAB  and  FORTRAN  IMSL. 

8.1  Background  IMSL  Optimizer 

The  IMSL  FORTRAN  library  contains  optimization  routines  classified  in  the  follow¬ 
ing  categories:  unconstrained  minimization,  minimization  with  simple  bounds,  linearly  con¬ 
strained  minimization,  and  nonlinearly  constrained  minimization.  The  routines  are  proprietary 
and  can  not  be  altered  by  the  user.  The  DNCONG  routine  in  that  library  solves  a  general 
nonlinear  programming  problem  using  the  SQP  algorithm  with  a  user  supplied  gradient.  It 
is  a  double  precision  routine.  The  optimization  algorithm  is  based  on  a  FORTRAN  code 
developed  by  Schittkowski  [Sch85].  First,  the  general  optimization  problem  is  stated. 


minimize 

subject  to: 

equality  constraints: 

gj(x)  =  0  i  =  l,...,me 

inequality  constraints: 

gj{x)  =  0  j  =  1 

side  constraints: 

Xi  <  X  <  Xu 

8.1.1  IMSL  Implementation  of  Quadratic  Sub-Problem .  The  IMSL  routine  includes 

two  improvements  that  make  the  code  more  efficient.  The  first  improvement  is  the  use  of 
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an  active  set  strategy  to  calculate  constraint  gradients.  This  reduces  the  number  of  gradient 
calculations  required.  The  active  set  flag  was  not  utilized  for  this  thesis;  therefore,  this  did  not 
effect  the  comparison  between  IMSL  and  MATLAB. 

The  second  improvement  in  IMSL’s  implementation  of  SQP  is  the  introduction  of  a 
variable  S  which  is  added  to  the  subproblem  to  minimize  possible  inconsistent  subproblems. 
This  equation  is  only  used  if  a  search  direction  cannot  be  found  for  ^  =  0.  The  QP  problem 
given  by  (2.104)  through  (2.106)  is  rewritten  here  as 


1  'T*  1  O 

minimize  s  V/  +  -s  VPs  -|-  -pk^ 

(8.1) 

subject  to 

/  _  \ 

(8.2) 

(8.3) 

\  ) 

gjixk)  +  s^Vgj{xkj)  >  0,  j  G  Kk 

(8.4) 

Xi  —  Xk  <  S  <  Xu  -  Xk 

(8.5) 

0<6<1 

(8.6) 

where  x  is  a  vector  of  design  variables  and  H  is  updated  using  a  quasi-Newton  method  to 
approach  the  Hessian  matrix.  The  update  method  used  in  IMSL  is  the  BFGS  method  discussed 
in  Chapter  2. 

8.1.2  IMSL  Implementation  of  the  Merit  Function.  The  merit  function  is  used  to 
determine  the  end  of  a  line  search.  The  next  values  of  the  design  variables  and  the  Lagrange 
multipliers  are  defined  at  the  end  of  the  line  search.  The  equations  for  the  design  variables 
(xfc+i)  and  Lagrange  multipliers  (vk+i)  are 


Xk+\  —  Xk  /3kSk 

(8.7) 

Vk+1  =Vk  +  ^k{uk  -  Vk) 

(8.8) 

where  Uk  is  the  Lagrange  multiplier  used  in  the  quadratic  sub-problem. 
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The  merit  function  is  minimized  so  that  the  design  vector  Xk+\  is  in  the  feasible  design 
space  and  the  design  vector  shows  some  improvement.  The  penalty  function  used  in  IMSL  is 


(8.9) 


where  is  a  penalty  function  and  is  a  vector  of  penalty  parameters.  The  IMSL  routine 
implements  the  augmented  Lagrangian  function  proposed  by  Schittkowski  as  the  penalty. 
This  method  has  been  determined  to  be  more  efficient  than  an  exact  penalty  function.  The 
augmented  Lagrange  penalty  is 


j=l  j=me+l 


-  0.5rjg]{x))  if  gj{x)  <  Vj/rj 

O.bVjjrj  otherwise 

(8.10) 


where  m'  =  m  +  2n  is  used  to  include  the  design  variable  bounds  as  constraints.  The  line 
search  starts  at  ^  =  1  and  is  reduced  until  (8.11)  is  met,  where  (^*,(0)<0  must  hold. 


Mot)  <  (pk{0)  + 


(8.11) 


where  ^  is  a  constant  and  (pk<0.  The  penalty  parameters  are  used  to  guarantee  convergence. 
There  is  a  different  penalty  parameter  for  each  constraint,  given  by 


=  max  ( 


—Vj  (fc))^ 
(l-5fc)s^tVs* 


),  j=l,...m 


(8.12) 


where 

<^r  =  min(l,  -^)  (8.13) 

V  ^ 

8.1 3  Convergence  Criteria.  The  convergence  criteria  for  IMSL  is  less  restrictive 
than  the  criteria  for  MATLAB.  It  is  not  an  option  to  input  to  the  IMSL  code.  The  source 
code  is  required  to  modify  this  parameter.  The  following  is  the  convergence  criteria  on  the 
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objective,  where  KTO  is  a  measure  of  the  Kuhn-Tucker  optimality  conditions. 


KTO  = 


j=i 


(8.14) 


The  convergence  criteria  on  constraints  is  measured  by  SCV  the  sum  of  constraint  values. 
This  criteria  is 

SCV  <  yfe  (8.15) 


The  variable  e  is  the  error  tolerance.  This  is  the  variable  that  the  user  can  set  if  the  source 
code  is  available. 


5.2  Background  MATLAB  Optimizer 

The  MATLAB  optimization  toolbox  contains  functions  that  solve  many  types  of  general 
nonlinear  optimization  problems.  Some  examples  of  optimization  problems  MATLAB  can 
solve  are  constrained  and  unconstrained  problems,  nonlinear  least  square  problems  and  matrix 
minimization  problems.  MATLAB  also  allows  the  user  to  supply  gradients.  If  gradients  are  not 
supplied  by  the  user,  finite  difference  gradients  are  calculated.  The  constrained  minimization 
routine  in  MATLAB  ’CONSTR’  uses  SQP.  The  MATLAB  implementation  is  discussed  in  this 
section. 

8.2.1  MATLAB  Implementation  of  Quadratic  Subproblem.  The  MATLAB  routine 
uses  a  quadratic  subproblem  to  calculate  the  search  direction.  It  is  a  quadratic  approximation 
of  the  Lagrangian  with  no  constant  term.  The  nonlinear  constraints  have  been  linearized  for 
this  problem.  The  quadratic  subproblem  that  the  MATLAB  routine  solves  is  given  in  (2.104) 
through  (2.106).  This  sub-problem  is  solved  using  a  Quadratic  Programming  (QP)  routine  in 
MATLAB.  The  solution  to  the  QP  problem  is  iterative.  The  MATLAB  QP  solver  calculates 
all  constraints,  decides  which  ones  are  active,  and  then  uses  those  to  solve  the  QP  problem. 
It  should  be  noted  that  aU  the  constraints  and  their  gradients  must  be  calculated,  and  then  the 
routine  determines  which  ones  are  active. 
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8.2.2  Merit  Function.  Once  the  search  direction  is  found  the  step  length  in  that 
direction  must  be  determined.  A  requirement  of  the  step  length  is  that  it  must  decrease  a  merit 
function  by  a  certain  amount.  The  MATLAB  implementation  has  two  merit  functions.  The 
first  merit  function  is  an  exact  penalty  that  is  written  as: 

me  m 

h  =  f{x)  +  m  i  0,  gi{x)  }  (8.16) 

1=1  i=me 


where 


Cki  maxAj,  '^1)5 i — l5...»m 


(8.17) 


This  penalty  function  allows  a  contribution  from  old  Lagrange  multipliers,  which  were  inactive 
in  the  QP  solution,  but  which  were  recently  active.  This  less  stringent  merit  function  was  added 
to  improve  convergence.  The  second  merit  function  is  used  to  improve  either  the  constraint  or 
the  objective  function  unless  the  QP  subproblem  is  infeasible.  If  the  QP  problem  is  infeasible 
the  routine  only  tries  to  reduce  the  maximum  constraint.  This  second  merit  function  is: 


maxj  gj{X)  if  X  is  infeasible 

if  X  infeasible  and  F(X)  >0  >  • 
F{X)  —  1  otherwise 


(8.18) 


Both  of  the  merit  functions  (8.16)  and  (8.18)  must  be  satisfied  and  the  number  of 
iterations  in  the  line  search  must  be  less  than  some  maximum  for  the  line  search  to  continue. 
If  a  step  length  of  one  does  not  meet  the  criteria  of  the  two  merit  functions,  the  step  length 
is  bisected  and  another  function  evaluation  is  done.  This  is  completed  up  to  the  number  of 
maximum  iterations.  Once  the  step  size  is  less  than  10"^,  the  step  length  is  set  negative  at 
every  other  bisection  step.  The  step  length  is  divided  fifteen  times  before  the  step  length  of 
10“^  is  set.  This  forces  the  optimizer  to  look  in  the  opposite  direction  than  the  original  search 
direction. 
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8.2.3  Convergence  Criteria.  The  convergence  criteria  for  the  objective,  constraint 
and  the  design  variables  for  MATLAB  can  be  input  by  the  user.  The  optimization  is  considered 
successful  when: 

1 .  The  maximum  absolute  change  in  the  search  direction  is  less  than  two  times  the  tolerance 
on  the  design  variables. 

2.  The  gradient  of  the  objective  multiplied  by  the  search  direction  is  less  than  two  times 
the  tolerance  on  the  objective. 

3.  The  maximum  constraint  is  less  than  or  equal  to  a  set  tolerance  from  zero  (i.e.,  the 
maximum  constraint  must  be  active).  This  was  modified  for  the  H2IH00  optimization 
problem  for  the  criteria  to  be  the  maximum  constraint  must  be  within  a  set  tolerance 
from  zero. 

83  Comparison  of  Methods 

The  IMSL  method  provides  an  improved  active  set  strategy  to  calculate  the  gradients 
compared  to  MATLAB.  This  means  that  fewer  gradients  should  be  calculated  by  the  user 
when  the  active  flag  is  used.  The  biggest  difference  in  this  comparison  is  that  IMSL  has  an 
improved  merit  function  and  line  search  methodology  compared  to  MATLAB.  Instead  of  a 
bisection  of  the  step  length  as  in  the  MATLAB  routine,  the  IMSL  interpolation  is  limited  by 
default  to,  at  most,  five  function  evaluations  per  line  search.  The  penalty  function  used  in 
the  merit  function  in  MATLAB  can  reduce  the  effectiveness  of  convergence  of  the  problem. 
Finally,  the  modified  QP  problem  in  IMSL  should  prevent  failure  to  produce  a  search  direction 
when  constraints  are  grossly  violated.  These  comparisons  lead  to  the  conclusions  that  IMSL 
should  be  more  efficient  and  more  accurate  in  the  optimization  process  than  the  MATLAB 
optimization  routine.  To  prove  this  conjecture,  the  F-16  SISO  example  was  run  using  both 
optimization  routines  and  comparisons  were  made.  This  is  discussed  in  Section  8.4. 
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Table  8 . 1  Convergence  Criteria  for  M ATLAB 


Comparison 

Objective 

Constraint 

1 

le-8 

le-4 

2 

le-8 

le-4 

3 

le-8 

le-4 

8.4  Comparison  Using  F-1 6  Example 

The  F-1 6  example  discussed  in  Chapter  4  was  run  using  the  IMSL/DNCONG  opti¬ 
mization  routine.  The  FORTRAN  program  called  MATLAB  as  an  "engine"  so  that  MATLAB 
routines  for  controls  were  available.  Therefore,  the  core  MATLAB  code  was  called  from  FOR¬ 
TRAN  to  generate  the  function  and  gradient  calculations  and  the  IMSL/DNCONG  routine 
was  used  for  the  H2IH00  optimization. 

Three  comparisons  were  made  for  this  example.  The  first  was  done  at  optimal  and 
then  reduced  in  7  increments  of  1000.  The  second  comparison  was  done  close  to  the  knee  of 
the  curve.  It  began  at  a  7  value  of  100  and  reduced  7  by  steps  of  5.  Finally,  a  comparison  was 
made  at  the  knee  of  the  mixed  curve  where  many  multiple  peaks  in  the  singular  value  plot 
were  evident.  This  comparison  started  at  a  7  level  of  5  and  was  reduced  by  7  steps  of  0.5. 
Note  the  absolute  7  for  a  seventh  order  controller  is  1 .32.  The  order  of  the  mixed  controller  for 
this  example  is  six.  These  three  comparisons  will  be  discussed  in  the  following  subsections. 

The  comparisons  were  made  as  close  as  possible  by  setting  the  tolerances  in  MATLAB 
as  close  to  the  tolerance  in  IMSL.  See  Tables  8.1  and  8.2  The  numbers  were  not  always 
exactly  the  same  since  the  criteria  for  IMSL  convergence  on  a  constraint  is  measuring  the  sum 
of  aU  the  violated  constraints  and  comparing  that  to  the  tolerance.  When  multiple  peaks  were 
constrained  and  caused  a  constraint  violation  this  tolerance  was  relaxed  to  provide  a  better 
comparison.  This  provided  as  close  a  match  for  convergence  between  the  two  routines  as 
possible.  The  maximum  number  of  iterations  within  a  line  search  was  also  adjusted  in  EMSL 
from  5  to  20. 
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Table  8.2  Convergence  Criteria  for  IMSL 


Comparison 

KTO 

scv 

1 

le-8 

le-4 

2 

le-8 

le-4 

3 

le-4 

le-2 

Figure  8.1  Singular  Value  Plot  MATLAB  (Comparison  1) 


8.4.1  Comparison  I:  Initial  7  at  H2  Optimal,  A'y=-1000  .  This  comparison 
was  done  to  investigate  the  performance  of  the  two  methods  beginning  at  the  H2  optimal 
compensator.  Large  steps  were  taken  to  identify  how  the  methods  handled  updating  the 
Hessian  matrix.  The  number  of  7  steps  was  set  at  70.  The  MATLAB  routine  became 
numerically  ill  conditioned  with  Hessian  update  errors  at  the  18th  7  step.  The  IMSL  routine 
completed  all  70  7  steps.  The  singular  value  plots  for  the  MATLAB  and  IMSL  examples  are 
shown  in  Figures  8.1  and  8.2. 

The  number  of  function  and  gradient  evaluations  is  presented  in  Table  8.3.  Only  the 
first  18  iterations  of  IMSL  are  presented  for  comparison  purposes,  since  that  was  the  number 
of  iterations  completed  by  MATLAB.  The  IMSL  routine  used  approximately  100  less  function 
and  gradient  evaluations  than  the  MATLAB  optimizer.  A  comparison  of  CPU  time  per  each 
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Figure  8.2  Singular  Value  Plot  IMSL  (Comparison  1) 


step  of  7  is  presented  in  Table  8.4,  where  the  IMSL  and  MATLAB  times  are  practically  the 
same.  This  was  true  even  though  the  number  of  function  and  gradient  evaluations  was  much 
less  for  IMSL.  This  may  be  due  to  more  time  spent  in  the  IMSL  optimizer  versus  time  spent  in 
the  MATLAB  optimizer  or  due  to  the  extra  overhead  for  IMSL  to  call  MATLAB  as  an  engine. 
Again,  the  first  18  steps  of  the  IMSL  run  are  compared  here  to  the  MATLAB  data. 

A  comparison  of  the  mixed  H2IH00  curves  are  given  in  Figure  8.3.  The  MATLAB 
curve  could  not  go  below  7=60,000.  It  is  difficult  to  see  on  the  plot,  but  the  MATLAB 
curve  at  7=61744  shows  a  dip  attributed  to  a  local  minimum  point.  This  dip  causes  the 
MATLAB  curve  to  be  located  on  the  IMSL  curve.  The  IMSL  mixed  curve  shows  a  smooth 
but  non-monotonically  increasing  plot. 

8.4.2  Comparison  2:  Initial  ')=100,  A'y=-5.  This  comparison  was  done  close  to 
the  knee  of  the  mixed  H2IH00  curve.  There  was  only  one  peak  for  most  of  the  iteration  so 
the  tolerance  used  was  10“®. 

Both  IMSL  and  MATLAB  completed  aU  19  iterations.  The  singular  value  plots  of 
MATLAB  and  IMSL  are  given  in  Figures  8.4  and  8.5.  The  MATLAB  optimizer  found  a 
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Table  8.3  Comparison  1  Number  of  Function  and  Gradient  Evaluation 


7 

MATLAB 

NUMFUN  NUMGRAD 

IMSL 

NUMFUN  NUMGRAD 

76744 

13 

12 

17 

15 

75744 

16 

12 

15 

10 

74744 

20 

14 

11 

9 

73744 

11 

9 

14 

8 

72744 

17 

9 

14 

9 

71744 

36 

19 

11 

7 

70744 

14 

6 

13 

8 

69744 

14 

6 

12 

7 

68744 

13 

6 

10 

6 

67744 

37 

19 

13 

7 

66744 

13 

6 

19 

14 

65744 

17 

11 

8 

4 

64744 

11 

5 

9 

5 

63744 

15 

9 

9 

5 

62744 

15 

10 

9 

5 

61744 

32 

22 

10 

5 

60744 

7 

5 

11 

6 

TOTAL  301  280  205  130 
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Table  8.4  Comparison  1  CPU  Time  of  MATLAB  vs  IMSL 


7 

MATLAB 

CPU  sec 

IMSL 

CPU  sec 

76744 

27.21 

31.92 

75744 

20.36 

21.87 

74744 

21.93 

19.00 

73744 

17.36 

19.33 

72744 

19.26 

21.07 

71744 

28.39 

18.47 

70744 

17.25 

20.55 

69744 

17.01 

19.70 

68744 

16.78 

18.36 

67744 

29.31 

21.20 

66744 

17.53 

27.54 

65744 

20.00 

16.34 

64744 

16.06 

17.64 

63744 

18.69 

19.11 

62744 

19.30 

18.55 

61744 

28.14 

19.00 

60744 

15.21 

20.22 

TOTAL 

349.88 

349.74 

different  design  than  the  IMSL  optimizer.  The  shapes  of  the  singular  value  plots  are  different 
at  the  first  values  of  7.  The  shape  of  the  MATLAB  singular  value  plot  indicates  a  local 
minimum  will  be  evident  on  the  mixed  H2IH00  curve  for  MATLAB  near  the  first  three  7 
steps.  The  local  minimum  for  the  IMSL  run  should  appear  at  7=30. 

The  number  of  function  and  gradient  evaluations  is  presented  in  Table  8.5.  The  number 
of  function  evaluations  for  IMSL  was  600  more  than  for  MATLAB.  This  increase  in  function 
evaluations  can  be  attributed  to  not  relaxing  the  tolerance  on  IMSL  enough  for  the  multiple 
peaks  to  match  the  tolerance  used  for  MATLAB.  Since  the  convergence  criteria  was  not  the 
same  for  MATLAB  and  IMSL  different  local  minima  were  found  during  the  optimization 
process.  This  is  evident  from  the  singular  value  plots  using  the  two  methods.  IMSL  did  not 
find  the  lower  local  minimum  as  soon  as  MATLAB  did.  Also,  due  to  different  convergence 
criteria,  the  shapes  of  the  singular  value  plots  are  different  between  7=90  through  7=30.  The 
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Singular  Value  Plot 


frequency 


Figure  8.4  Singular  Value  Plot  MATLAB  (Comparison  2) 

singular  Value  Plot 


Figure  8.5  Singular  Value  Plot  IMSL  (Comparison  2) 


time  in  CPU  for  both  methods  is  shown  in  Table  8.6.  The  increase  in  time  is  due  to  the  increase 
in  the  number  of  function  evaluations.  This  increase  can  also  be  attributed  to  the  tolerance  set 
for  the  convergence  criteria  of  the  IMSL  routine. 

The  mixed  H2  /  Hoo  curves  show  a  significant  difference  in  local  minima.  The  MATLAB 
optimization  routine  found  a  local  minimum  after  three  7  steps.  The  IMSL  routine  found  a 
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Table  8.5  Comparison  2  Number  of  Function  and  Gradient  Evaluation 


MATLAB  IMSL 


7 

NUMFUN 

NUMGRAD 

NUMFUN 

NUMGRAD 

95 

41 

20 

190 

100 

90 

28 

10 

13 

8 

85 

178 

118 

9 

7 

80 

12 

8 

10 

8 

75 

9 

4 

9 

7 

70 

10 

6 

9 

7 

65 

6 

6 

12 

10 

60 

9 

4 

23 

15 

55 

8 

6 

20 

12 

50 

7 

4 

34 

23 

45 

7 

6 

17 

11 

40 

11 

7 

12 

8 

35 

5 

4 

14 

12 

30 

4 

4 

186 

100 

25 

4 

4 

147 

81 

20 

5 

5 

5 

5 

15 

4 

4 

7 

7 

10 

8 

8 

12 

10 

5 

8 

8 

268 

98 

TOTAL 

364 

236 

997 

529 

local  minimum  after  13  steps  in  7.  This  difference  in  local  minimum  can  be  attributed 
to  the  difference  in  convergence  criteria  used  between  the  two  methods.  Both  curves  do 
end  at  the  same  mixed  design  points  but  both  took  different  paths.  Also  when  comparing 
function  evaluations  and  time  required  it  can  be  seen  that  after  MATLAB  found  the  lower 
local  minimum  it  was  much  quicker  in  calculating  the  rest  of  the  curve.  IMSL,  on  the  other 
hand,  took  longer  to  find  that  local  minimum  and  required  more  function  evaluations  and  time. 
The  MATLAB  and  IMSL  mixed  curves  are  presented  in  Figure  8.8. 

Since  it  is  evident  that  the  IMSL  routine  was  not  tracking  the  same  local  minimum 
as  MATLAB  due  to  convergence  criteria  the  comparison  is  not  valid.  It  would  be  a  better 
comparison  for  number  of  function  evaluations  and  time  if  both  MATLAB  and  IMSL  were 
tracking  the  same  local  minimum.  Removing  the  scale  factor  of  (0.0001)  from  the  objective 
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Alpha 


Table  8.6  Comparison  2  CPU  Time  of  MATLAB  vs  IMSL 


7 

MATLAB 
CPU  sec 

IMSL 
CPU  sec 

95 

43.95 

153.33 

90 

27.41 

20.05 

85 

148.86 

17.93 

80 

20.53 

18.56 

75 

16.84 

18.11 

70 

18.66 

18.01 

65 

17.04 

20.43 

60 

17.45 

29.38 

55 

18.18 

26.65 

50 

16.16 

37.21 

45 

17.05 

25.13 

40 

19.63 

20.03 

35 

17.66 

24.18 

30 

16.26 

166.05 

25 

16.69 

249.13 

20 

18.03 

18.15 

15 

16.45 

21.03 

10 

22.59 

24.86 

5 

22.81 

260.88 

TOTAL 

512.35 

1169.1 

Figure  8.6  Mixed  H2IH00  Curve  MATLAB-’o’  IMSL-’x’  (Comparison  2) 
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Singular  Value  Plot 


trequency 


Figure  8.7  Singular  Value  Plot  No  Scale  IMSL  (Comparison  2) 

function  for  IMSL,  caused  IMSL  to  find  the  same  local  minimum  as  MATLAB  found  much 
earlier.  Since  the  scale  factor  on  the  objective  was  removed,  the  convergence  criteria  was  set 
to  10“^  to  allow  IMSL  to  converge  in  20  steps  in  the  line  search.  A  plot  of  the  singular  value 
plot  using  IMSL  with  no  scale  factor  on  the  objective  is  shown  in  Figure  8.7.  IMSL  then  took 
one  7  step  to  achieve  the  same  singular  value  plot  shape  as  MATLAB.  A  plot  of  the  mixed 
H2  / Hoo  curve  with  the  MATLAB  curve  and  the  IMSL  curve  with  no  scale  on  the  objective  is 
shown  in  Figure .  IMSL  and  MATLAB  compare  much  better.  Finally,  without  the  scale  factor 
for  IMSL  tables  of  function  evaluations  and  time  comparisons  are  shown  in  8.7  and  8.8. 

Notice  that  even  though  the  scale  factor  has  been  taken  off  the  objective  for  IMSL,  the 
IMSL  runs  took  less  time  and  less  function  evaluations  than  the  previous  IMSL  run  with  the 
scales.  This  can  be  attributed  to  finding  the  lower  local  minimum  earlier  with  the  scale  than 
without  the  scale.  The  final  conclusion  from  this  comparison  is  that  without  the  scale  factor 
on  the  objective  for  IMSL,  the  mixed  curve  for  MATLAB  and  IMSL  compared  very  well.  The 
time  and  number  of  function  evaluations  is  still  more  for  IMSL  than  MATLAB. 

8.43  Comparison  3:  Initial  Gamma=5,  A‘y=-0.25.  The  final  comparison  is  very 
close  to  Hoo  optimal  for  this  problem.  Again,  7  for  this  problem  is  1.32.  This  comparison  is 
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1 

MATLAB 

NUMFUN  NUMGRAD 

IMSL 

NUMFUN  NUMGRAD 

95 

41 

20 

174 

88 

90 

28 

10 

27 

11 

85 

178 

118 

103 

28 

80 

12 

8 

29 

12 

75 

9 

4 

27 

10 

70 

10 

6 

24 

12 

65 

6 

6 

26 

10 

60 

9 

4 

27 

12 

55 

8 

6 

28 

11 

50 

7 

4 

37 

20 

45 

7 

6 

27 

13 

40 

11 

7 

25 

12 

35 

5 

4 

27 

15 

30 

4 

4 

37 

19 

25 

4 

4 

41 

20 

20 

5 

5 

41 

25 

15 

4 

4 

28 

16 

10 

8 

8 

27 

16 

5 

8 

8 

33 

24 

TOTAL 

364 

236 

840 

374 
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Table  8.8  Comparison  2  CPU  Time  of  MATLAB  vs  IMSL 


7 

MATLAB 

CPU  sec 

IMSL 

CPU  sec 

95 

43.95 

153.43 

90 

27.41 

30.20 

85 

148.86 

86.90 

80 

20.53 

32.38 

75 

16.84 

30.96 

70 

18.66 

29.76 

65 

17.04 

30.16 

60 

17.45 

32.83 

55 

18.18 

31.83 

50 

16.16 

42.68 

45 

17.05 

34.03 

40 

19.63 

32.05 

35 

17.66 

40.76 

30 

16.26 

50.26 

25 

16.69 

53.30 

20 

18.03 

57.65 

15 

16.45 

42.20 

10 

22.59 

42.30 

5 

22.81 

49.10 

TOTAL 

512.35 

909.7 

dov^^n  to  7  =  1.5.  Due  to  four  peaks  being  constrained  at  the  start  of  the  example  and  up  to  10 
peaks  total  per  iteration  being  constrained,  the  IMSL  convergence  tolerance  was  relaxed  to 
10“^.  This  was  the  best  comparison  possible  for  convergence  criteria.  Both  cases  completed 
aU  19  iterations.  The  singular  value  plots  for  MATLAB  and  IMSL  are  approximately  the  same 
shape  except  at  the  frequencies  between  10“®  and  10“^.  The  IMSL  plot  shows  a  different 
design  being  caused  by  a  local  minimum  at  the  lowest  value  of  7.  However,  the  lowest  7 
singular  value  curve  matches  for  IMSL  and  MATLAB.  They  are  shown  in  Figures  8.9  and 
8.10. 

The  number  of  function  and  gradient  evaluations  is  presented  in  Table  8.9.  The  IMSL 
routine  required  200  less  function  evaluations  and  100  less  gradient  evaluations  than  MAT- 
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6, 


Singular  Value  Plot 


Figure  8.9  Singular  Value  Plot  MATLAB  (Comparison  3) 


Figure  8.10  Singular  Value  Plot  IMSL  (Comparison  3) 


LAB.  A  comparison  of  CPU  time  for  each  step  of  7  is  presented  in  Table  8.10.  The  time  was 
reduced  by  over  1000  CPU  seconds  by  using  IMSL. 

The  mixed  H^jHoo  curves  are  shown  in  Figure  8.11.  The  local  minimum  dip  in  the 
mixed  curve  for  IMSL  is  evident  at  the  lowest  7  level  found.  The  two-norm  and  infinity-norm 
match  at  the  lowest  7  value  for  both  the  mixed  curves  of  IMSL  and  MATLAB. 


Table  8.9  Comparison  3  Number  of  Function  and  Gradient  Evaluation  Comparison 


MATLAB  IMSL 


7 

NUMFUN 

NUMGRAD 

NUMFUN 

NUMGRAD 

4.75 

18 

18 

2 

2 

4.50 

6 

6 

2 

2 

4.25 

4 

4 

2 

2 

4.00 

17 

17 

2 

2 

3.75 

6 

6 

5 

4 

3.50 

6 

6 

3 

3 

3.25 

7 

7 

10 

6 

3.00 

37 

32 

6 

5 

2.75 

22 

15 

100 

35 

2.50 

32 

26 

25 

11 

2.25 

32 

27 

28 

12 

2.00 

204 

37 

38 

16 

1.75 

155 

64 

55 

21 

1.5 

42 

21 

63 

26 

TOTAL 

588 

286 

341 

147 

Table  8.10  Comparison  3  CPU  Time  of  MATLAB  vs  IMSL 

7 

MATLAB 

CPU  sec 

IMSL 
CPU  sec 

4.75 

49.64 

19.78 

4.50 

19.33 

13.22 

4.25 

16.48 

13.06 

4.00 

34.73 

13.44 

3.75 

19.85 

17.60 

3.50 

19.64 

15.17 

3.25 

20.90 

23.96 

3.00 

71.76 

20.21 

2.75 

44.44 

133.05 

2.50 

63.65 

50.43 

2.25 

73.30 

56.25 

2.00 

305.15 

81.39 

1.75 

422.43 

102.31 

1.5 

465.20 

127.57 

TOTAL 

1626.50 

687.44 
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Figure  8.11  Mixed  H2I Hoc,  Curve  MATLAB-’o’  IMSL-’x’  (Comparison  3) 

8.5  Summary  and  Conclusions 

The  number  of  function  and  gradient  evaluations  for  each  of  the  three  cases  is  signif¬ 
icantly  less  by  approximately  (100  to  200  evaluations)  for  IMSL  than  MATLAB.  The  time 
required  to  calculate  the  mixed  controller  is,  in  general,  less  for  IMSL  than  MATLAB.  This 
can  be  attributed  to  less  function  and  gradient  calculations  and  the  use  of  FORTRAN  code 
which  is,  in  general,  faster  than  MATLAB.  The  singular  value  plots  were  shown  to  compare 
the  general  trends  of  the  peaks  using  each  method.  They  were  generally  the  same  using  either 
method,  but  there  were  examples  of  different  designs  being  found  by  either  optimizer.  The 
mixed  H2  /  ifoo  curves  were  generally  smoother  for  IMSL. 

The  first  comparison  showed  the  IMSL  curve  to  be  smooth  with  no  local  minimum 
jumps.  However,  it  was  not  monotonically  increasing.  The  MATLAB  curve  showed  a 
monotonically  increasing  curve  with  a  local  minimum  dip  at  the  left  of  the  curve.  The  values 
of  the  two-norm  were  approximately  the  same  in  the  range  of  7  that  could  be  compared. 
The  second  comparison  showed  IMSL’s  mixed  curve  as  a  smooth  monotonically  increasing 
function.  The  MATLAB  plot  showed  significant  local  minimum  dips  in  the  curve.  The  other 
observation  is  that  the  MATLAB  curve  demonstrated  lower  values  of  the  two-norm  for  7 
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between  (5  to  180).  The  last  comparison  demonstrated  very  similar  mixed  H2IH00  curves  for 
both  MATLAB  and  IMSL.  They  were  both  monotonically  increasing  with  approximately  the 
same  two-norm  value  at  different  7  levels.  The  IMSL  plot  did  show  a  dip  in  the  mixed  curve 
due  to  a  local  minimum,  but  ended  at  the  same  design  point  as  MATLAB. 

The  conclusion  from  this  comparison  is  that  the  IMSL/DNCONG  routine  is  a  better 
optimizer  than  the  MATLAB  ’CONSTR’  function  for  mixed  H2IH00  optimization.  It  is  faster 
and  takes  less  function  and  gradient  evaluations.  It  also  helps  to  calculate  a  smoother  mixed 
H2IH00  curve  but  does  demonstrate  local  minimum  jumps.  The  second  comparison,  where 
the  MATLAB  curve  showed  a  lower  two-norm,  means  that  the  MATLAB  program  found  a 
lower  local  minimum  than  the  IMSL  program.  Since  this  is  not  a  convex  problem,  this  is  to 
be  expected. 

The  final  conclusion  from  this  chapter  is  that  the  mixed  H2/H00  problem  should  be 
run  using  the  Shittkowski  penalty  method  for  the  merit  function  and  his  method  to  update  the 
quadratic  sub-problem.  These  are  implemented  in  IMSL  but  are  not  easily  accessible  to  the 
user  to  make  modifications  to  tolerances  and  number  of  evaluations.  Also,  calling  MATLAB 
from  FORTRAN  is  cumbersome  to  code  and  very  difficult  to  debug.  The  best  solution  would 
be  to  implement  the  Shittkowski  methods  in  place  of  MATLAB ’s  ’CONSTR’  function.  This 
would  improve  the  speed  and  convergence  of  the  mixed  H2IH00  problem. 
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IX.  Conclusions  and  Recommendations 


9.1  Thesis  Summary  and  Conclusions 

This  thesis  identified  ways  to  improve  the  run  time  and  the  reliability  of  the  mixed 
H2IH00  optimization  process.  The  improved  code  was  run  for  a  realistic  flight  control 
problem  of  an  F-16.  The  controller  design  for  the  realistic  F-16  problem  was  done  at  fifth, 
sixth  and  seventh  order.  Comparisons  were  then  made  between  the  mixed  H2IH00  controllers 
at  different  orders  to  the  LQG/LTR  design.  The  improved  code  was  then  used  to  generate 
mixed  H2/H00  controllers  for  a  MIMO  missile  problem.  Mixed  H2IH00  curves  were  found 
for  both  the  3x3  MIMO  and  the  5x3  MIMO  designs.  Finally,  a  comparison  of  the  MATLAB 
’CONSTR’  optimization  function  and  the  FORTRAN  IMSL/DNCONG  optimization  routine 
was  done. 

The  purpose  of  this  thesis  was  to  improve  the  run  time  and  reliability  of  the  mixed 
H2  /  Hoo  optimization  process  and  then  demonstrate  its  effectiveness  on  realistic  control  prob¬ 
lems.  The  examples  discussed  show  that  the  run  time  and  reliability  of  generating  a  mixed 
H2I Hoo  controller  has  significantly  improved.  The  results  from  this  thesis  show  that  the 
mixed  H2I  Hoo  design  is  better  for  the  F-16  example  than  the  LQG/LTR  example  even  for 
a  controller  that  is  three  orders  less  than  the  LQG/LTR.  It  has  also  demonstrated  that  where 
numeric  problems  for  the  MIMO  missile  example  prevented  mixed  H2/ Hoo  designs,  these 
controllers  can  now  be  found.  Finally,  the  comparison  between  the  MATLAB  and  IMSL  op¬ 
timization  routines  proved  that  the  IMSL  optimization  method  was  better  for  mixed  H2I Hoo 
optimization. 

9.2  Recommendations 

This  study  has  raised  several  new  questions  about  the  mixed  H2I Hoo  optimization 
design  process.  The  following  are  suggested  avenues  of  study  which  the  results  of  this  thesis 
indicate  to  be  worthwhile. 
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•  Code  Shittkowski’s  optimization  methods  in  MATIAB.  This  would  provide  the  designer 
with  the  efficiency  of  the  IMSL  optimization  routine  and  the  user  friendliness  of  MAT- 
LAB.  It  would  also  provide  the  ability  to  make  modifications  to  the  optimization  code 
when  necessary,  such  as  the  change  to  only  evaluate  the  maximum  Hoc  constraint. 

•  Investigate  way  to  reduce  the  probability  of  finding  a  local  minimum  in  the  mixed 
H2IH00  optimization  process.  There  are  methods  called  simulated  annealing  and 
genetic  algorithms  that  could  be  added  to  this  code  to  improve  the  chances  of  the 
optimization  routine  moving  away  from  a  local  minimum.  Implementing  these  meth¬ 
ods  could  improve  the  mixed  H2/H00  curve  by  improving  the  chances  that  a  global 
minimum  could  be  found. 

•  Scale  the  design  variables.  This  would  be  done  to  possibly  increase  efficiency.  The 
functions,  constraints  and  their  gradients  were  all  scaled  for  this  study.  The  other 
possibility  to  improve  performance  of  the  optimization  process  is  to  scale  the  design 
variables. 
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Appendix  A.  SISO  Example:  Matrices  for  Underlying  H2  and  Hqo  Problems 


The  matrices  for  the  SISO  example’s  underlying  H2  problem  are: 


-1.4850E-  02 

3.7382E+  01 

-3.2200E  +  01 

-1.7940E  +  01 

2.1400E-03 

O.OOOOE  +  00 

-8.0000E  -  05 

-1.4910E  +  00 

-1.3000E-03 

9.9600E-  01 

-1.8800E-  01 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

l.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-3.6000E  -  04 

9.7530E+  00 

2.9000E  -  04 

-9.6000E  -  01 

-1.9040E  +  01 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

-2.0000E  +  01 

O.OOOOE  +  00 

1.5000E-  03 

3.5264E+  01 

2.7200E-02 

-3.3400E  -  01 

-4.3660E  +  00 

-4.0000E  +  01 

8.3589E  -  01 

-3.3340E  -  02 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

Byj  — 

2.1808E  -  01 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE +00 

B 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
2.0000E  +  01 
O.OOOOE  +  00 


-1.5000E-  03 
O.OOOOE +00 


-3.5264E  +  01 
O.OOOOE  +  00 


-2.7200E  -  02 
O.OOOOE +00 


3.3400E  -  01 
O.OOOOE +00 


4.3660E  +  00 
O.OOOOE  +  00 


8.0000E+  01 
O.OOOOE +00 


Dzw 


0.0000E  +  00  O.OOOOE  +  00 
O.OOOOE  +  00  O.OOOOE  +  00 
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Dzu 


O.OOOOE  +  00 
l.OOOOB+01 


C 


y2 


-1.5000E-03  -3.5264E  +  01  -2.7200E-02  3.3400E-  01  4.3660E  +  00  8.0000E+01 


Dyui 


O.OOOOE  +  00  4.0000E  -  03 


Dyu  -  O.OOOOE +00 


The  matrices  for  the  SISO  example’s  underlying  Hoo  problem  are: 


-1.4850E-02 

3.7382E+01 

-3.2200E+01 

-1.7940E+01 

2.1400E-  03 

O.OOOOE +  00 

O.OOOOE  + 

-8.0000E  -  05 

-1.4910E+00 

-1.3000E-03 

9.9600E  -  01 

-1.8800E-  01 

O.OOOOE +00 

O.OOOOE  + 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

1.0000E+  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  + 

-3.6000E  -  04 

9.7530E+  00 

2.9000E-04 

-9.6000E  -  01 

-1.9040E  +  01 

O.OOOOE +00 

O.OOOOE  + 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

-2.0000E  +  01 

O.OOOOE +  00 

O.OOOOE  + 

1.5000E-03 

3.5264B+01 

2.7200E-02 

-3.3400E  -  01 

-4.3660E  +  00 

-4.0000E  +  01 

O.OOOOE  + 

-1.5000E-  03 

-3.5264E+01 

-2.7200E-02 

3.3400E-  01 

4.3660E  +  00 

8.0000E+01 

-l.OOOOE- 

J5d 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +00 
O.OOOOE  +  00 
l.OOOOE+00 
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O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +00 
=  O.OOOOE +00 
2.0000E  +  01 
O.OOOOE  +  00 
O.OOOOE +  00 


-1.5000E-  03  -3.5264E+01  -2.7200E-02  3.3400E  -  01  4.3660E+00  8.0000E  +  01  l.OOOOE+00 


Ded  -  I  l.OOOOE+00 

Ceu  =  O.OOOOE +00 

-1.5000E-03  -3.5264E+01  -2.7200E-02  3.3400E-01  4.3660E+00  8.0000E+01  O.OOOOE +00 


l.OOOOE+00 


Appendix  B.  3x3  MIMO  Problem:  Matrices  for  Underlying  H2  and  Hqo 

Problems 

The  matrices  for  the  3x3  MIMO  example’s  underlying  Hi  problem  are: 


-4.0000E  +  02 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

-4.0000E  +  02 

O.OOOOE+OO 

O.OOOOE+00 

O.OOOOE+OO 

-4.0000E  +  02 

-8.7451E-01 

1.3044E+02 

-4.5412E  -  01 

-2.4756E-01 

6.2620E  -  01 

-1.2753E  +  02 

4.5574E+06 

1.2009E  +  06 

-4.2461E  +  03 

-1.0943E  +  03 

1.2844E  +  05 

5.8025E  +  01 

-1.2180E  +  02 

5.6387E+02 

1.2598E  +  05 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  OO 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  + 

00 

-9.8435E  -  01 

-9.2340E-  02 

O.OOOOE+OO 

l.OOOOE+00 

O.OOOOE  + 

00 

-9.2340E  -  02 

-9.8435E-01 

O.OOOOE  +  00 

O.OOOOE+OO 

-1.0000E  + 

00 

-2.6737E+02 

2.6737E+  02 

-1.5949E+  00 

O.OOOOE+OO 

O.OOOOE  + 

00 

-1.9462E+  02 

-6.9672E  -  01 

O.OOOOE  +  00 

-1.5498E+00 

O.OOOOE  + 

00 

6.9673E  -  01 

1.9462E+02 

O.OOOOE  +  00 

O.OOOOE+OO 

-1.5498E  + 

00 

O.OOOOE  +  OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

-1.9687E-  03 

-1.8468E-05 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  + 

00 

-1.8468E-04 

-1.9687E-04 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  + 

00 

-5.3474E  -  01 

5.3474E  -  02 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

-3.8924E  -  01 

-1.3934E-04 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  + 

00 

1.3935E-  03 

3.8924E  -  02 

O.OOOOE+OO 

O.OOOOE  +  OO 

O.OOOOE  + 

00 
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l.OOOOE  +  OO 

O.OOOOE  +  00 

O.OOOOE  +  OO 

O.OOOOE  +  00 

l.OOOOE  +  OO 

O.OOOOE+OO 

O.OOOOE+00 

O.OOOOE  +  00 

l.OOOOE+OO 

^ti2  — 

O.OOOOE+00 

O.OOOOE  +  00 

O.OOOOE  +  OO 

O.OOOOE  +  OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

-S.4953E  +  02 

1.3982E+03 

-2.8476E  +  05 

-1.9504E  +  03 

2.9126E+05 

-1.0140E+  03 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  OO 

l.OOOOE  +  OO 

O.OOOOE+OO 

O.OOOOE  + 

00 

-2.0619E+02 

-2.1979E  +  03 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  + 

00 

-2.1979E+  03 

-2.0619E  +  02 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  OO 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  +  OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  + 

00 

_  O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

Dzu  — 

O.OOOOE  +  00 

7.071  IE +01 

O.OOOOE  +  00 

O.OOOOE  +  OO 

O.OOOOE+OO 

O.OOOOE  +  OO 

O.OOOOE  +  00 

7.071  IE +  01 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

7.071  IE +01 
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Cy,  = 

O.OOOOE +00 

-5.4953E  +  02 

-1.9504E  +  03 

O.OOOOE  +  00 

1.3982E  +  03 

2.9126E  +  05 

O.OOOOE  +  00 

-2.8476E+  05 

-1.0140E+03 

O.OOOOE +00 

-2.0619E+02 

-2.1979E+03 

O.OOOOE +00 

-2.1979E  +  03 

-2.0619E  +  02 

l.OOOOE+00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +  00 

Dyw  = 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

5.4538E  -  05 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

1.2500E-03 

O.OOOOE  +  00 

Dyu  — 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE +00 
O.OOOOE +00 
1.2500E-  03 
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The  matrices  for  the  3x3  MIMO  example’s  underlying  Hoc  problem  are 


-4.0000E  +  02 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE  +  00 

-4.0000E+  02 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-4.0000E+02 

O.OOOOE +00 

O.OOOOE +  00 

-8.7451E  -  01 

1.3044B+02 

-4.5412E-  01 

-9.8435E  -  01 

-9.234E  -  02 

-2.4756E-01 

6.2620E  -  01 

-1.2753E+02 

-9.2340E  -  02 

-9.8435E  -  01 

4.5574E  +  06 

1.2009E+  06 

-4.2461B+03 

-2.6737E  +  02 

2.6737E+02 

-1.0943E  +  03 

1.2844E+  05 

5.8025E+01 

-1.9462E  +  02 

-6.9672E  -  01 

-1.2180E  +  02 

5.6387E  +  02 

1.2598E+05 

6.0673E  -  01 

1.9462E+02 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +  00 

-5.4736E  +  02 

1.3982E+03 

-2.8476E+05 

-2.0619E+  02 

-2.1979E  +  03 

-1.9502B  +  03 

2.9126E+  05 

-1.0140E+03 

-2.1979E+  03 

-2.0619E  +  02 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

l.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

-l.OOOOE  +  OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-1.5949E+00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

-1.5498E  +  00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

-1.5498E  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

l.OOOOE+00 

O.OOOOE +  00 

O.OOOOE +00 

-l.OOOOE-  04 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-l.OOOOE  -  04 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-l.OOOOE  -  04 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

l.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

l.OOOOE+OO 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

l.OOOOE  +  OO 

B-4 


B 


«oo 


1.0000E+  00 
O.OOOOE+00 
O.OOOOE+00 
O.OOOOE+OO 
O.OOOOE+00 
O.OOOOE  +  00 
I5.0000E+00 
O.OOOOE+OO 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE+OO 


O.OOOOE  +  00 
l.OOOOE  +  OO 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE+OO 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE  +  00 
O.OOOOE+OO 
l.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  + 

Ce  = 

-5.4953E  + 

02 

1.3982E+03 

-2.8476E+05 

-2.0619E+  02 

-2.1979E  + 

-1.9502E  + 

03 

2.9126E  +  05 

-1.0140E+03 

-2.1979E+03 

-2.0619E  + 

l.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE  +  00 


O.OOOOE  +  00 
O.OOOOE+OO 
O.OOOOE  +  00 


O.OOOOE+OO 
O.OOOOE  +  OO 
O.OOOOE+OO 


2.0000E  +  01 
O.OOOOE+OO 
O.OOOOE+OO 


O.OOOOE  +  00 
2.0000E  +  01 
O.OOOOE  +  00 


O.OOOOE  +  00 
O.OOOOE+OO 
2.0000E+01 


C 


y<x 


l.OOOOE  +  OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

l.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

l.OOOOE  +  OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  + 

00 

-5.4736E+02 

1.3982E+03 

-2.8476E  +  05 

-2.0619E  +  02 

-2.1979E  + 

03 

-1.9502E  +  03 

2.9126E+05 

-1.0140E  +  03 

-2.1979E  +  03 

-2.0619E  + 

02 

l.OOOOE  +  OO 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE  +  OO 
O.OOOOE  +  00 
O.OOOOE  +  OO 


O.OOOOE  +  00 
O.OOOOE+OO 
O.OOOOE  +  00 


O.OOOOE+OO 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE+OO 
O.OOOOE  +  00 
O.OOOOE+OO 


O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 


DyJl 


l.OOOOE+OO 
O.OOOOE  +  00 
O.OOOOE+OO 


O.OOOOE  +  00 
l.OOOOE  +  OO 
O.OOOOE  +  00 


O.OOOOE+OO 
O.OOOOE  +  00 
l.OOOOE+OO 


B-5 


B-6 


Appendix  C.  5x3  MIMO  Problem:  Matrices  for  Underlying  H2  and 

Problems 

The  matrices  for  the  5x3  MIMO  example’s  underlying  H2  problem  are: 


-4.0000E+  02 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00  .. 

O.OOOOE  +  00 

-4.0000E  +  02 

O.OOOOE  +  00 

O.OOOOE +00  .. 

O.OOOOE +00 

O.OOOOE +00 

-4.0000E  +  02 

O.OOOOE +00  .. 

^2  = 

-8.7349E  -  01 

-2.4611E-  01 

1.3044E+02 

6.2620E-01 

-4.5412E  -  01 

-1.2753E  +  02 

-9.8435E  -  01  . . 

-9.2340E  -  02  . . 

3.3304E+06 

1.2009E+06 

-4.2461E  +  03 

-2.6737E  +  02  .. 

-1.0950E+  03 

1.2844E+05 

5.8025E  +  01 

-1.9462E  +  02  .. 

-1.2044E+  02 

5.6387E  +  02 

1.2598E+05 

6.9673E-01  .. 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-9.2340E  -  02 

O.OOOOE +  00 

l.OOOOE+00 

O.OOOOE  +  00 

-9.843SE  -  01 

O.OOOOE +00 

O.OOOOE  +  00 

-l.OOOOE+OO 

2.6737E  +  02 

-1.5949E+00 

O.OOOOE +00 

O.OOOOE +00 

-6.9672E  -  01 

O.OOOOE +00 

-1.5498E+00 

O.OOOOE +00 

1.9462E  +  02 

O.OOOOE  +  00 

O.OOOOE +00 

-1.5498E+  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  + 

00 

-1.9687E-  03 

-1.8468E-05 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

-1.8468E-  04 

-1.9687E-04 

O.OOOOE +  00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  + 

00 

-5.3474E-  01 

5.3474E-02 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  + 

00 

-3.8924E  -  01 

-1.3934E  -  04 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  + 

00 

1.3935E-  03 

3.8924E-02 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  + 

00 

l.OOOOE+OO 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

l.OOOOE+OO 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

l.OOOOE+OO 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

— 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-5.4953E  +  02 

1.3982E+  03 

-2.8476E+  05 

-1.9504E  +  03 

2.9126E+05 

-1.01 40E+  03 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

l.OOOOE  +  OO 

O.OOOOE +00 

O.OOOOE  + 

00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

l.OOOOE  +  OO 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE +  00 

1.0000E  + 

00 

-2.0619E+  02 

-2.1979E+03 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE  + 

00 

-2.1979E+  03 

-2.0619E+02 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE  + 

00 

O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +00 
O.OOOOE  +  00 
O.OOOOE +00 


O.OOOOE  +  00 
O.OOOOE +  00 
O.OOOOE +00 
O.OOOOE  +  00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +  00 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +  00 
O.OOOOE +  00 
O.OOOOE  +  00 


O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE  +  00 
O.OOOOE +  00 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +  00 
O.OOOOE  +  00 
O.OOOOE +  00 
O.OOOOE +  00 


C-2 


Dzu 


O.OOOOE  +  00 
O.OOOOE+00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
7.071  IE +  01 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
7.0711E  +  01 
O.OOOOE  +  00 


O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
7.0711E+01 


C 


Vi 


O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

-5.4953E  +  02 

1.3982E+03 

-2.8476E+05 

-1.9504E  +  03 

2.9126E+05 

-1.0140E+03 

O.OOOOE  +  00 

O.OOOOE  +  00 

l.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

l.OOOOE+00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

l.OOOOE+00 

-2.0619E+02 

-2.1979E+03 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

-2.1979E+03 

-2.0619E+02 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE  +  00 
O.OOOOE  +  00 
Oyiu  =  O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE +00 


O.OOOOE  +  00  5.4538E  -  05 

O.OOOOE  +  00  O.OOOOE  +  00 
O.OOOOE  +  00  O.OOOOE  +  00 
O.OOOOE  +  00  O.OOOOE  +  00 
O.OOOOE  +  00  O.OOOOE  +  00 


O.OOOOE  +  00  O.OOOOE  +  00 
5.4538E  -  05  O.OOOOE  +  00 
O.OOOOE  +  00  5.4538E  -  05 

O.OOOOE  +  00  O.OOOOE  +  00 
O.OOOOE  +  00  O.OOOOE  +  00 


O.OOOOE  +  00  O.OOOOE  +  00 
O.OOOOE  +  00  O.OOOOE  +  00 
O.OOOOE  +  00  O.OOOOE  +  00 
4 .0250E  -  02  O.OOOOE  +  00 
O.OOOOE  +00  4 .0250E  -  02 


Dyu 


O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 
O.OOOOE  +  00 


O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE +00 
O.OOOOE  +  00 


O.OOOOE  +  00 
O.OOOOE +00 
O.OOOOE  +  00 
O.OOOOE +00 
O.OOOOE +00 


C-3 


The  matrices  for  the  3x5  MEMO  example’s  underlying  Hoo  problem  are: 


-4.0000E  + 

02 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE  + 

00 

-4.0000E+02 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE +  00 

0.0000E  + 

00 

O.OOOOE +00 

-4.0000E+02 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

-8.7349E  - 

01 

1.3044E+  02 

-4.5412E-  01 

-9.8435E  -  01 

-9.2340E  -  02 

O.OOOOE +00 

-2.4611E- 

01 

6.2620E-  01 

-1.2753E+02 

-9.2340E  -  02 

-9.8435E  -  01 

O.OOOOE +00 

3.3304E  4- 

06 

1.2009E+  06 

-4.2461E+03 

-2.6737E+  02 

2.6737E  +  02 

-1.5949E  +  00 

-1.0950E  + 

03 

1.2844E+05 

5.8025E+01 

-1.9462E+  02 

-6.9672E  -  01 

O.OOOOE +00 

-1.2044E  + 

02 

5.6387E  +  02 

1.2598E+05 

6.9673E  -  01 

1.9462E+02 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

2.0000E  +  01 

O.OOOOE  + 

00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

-1.0991E  + 

04 

2.7964E+04 

-5.6952E+06 

-4.1238E+  03 

-4.3958E  +  04 

O.OOOOE +00 

-3.9008E  + 

04 

5.8252E+06 

-2.0280E+04 

-4.3958E  +  04 

-4.1238E  +  03 

O.OOOOE +  00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

1.0000E  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE +00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

-l.OOOOE+00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

-1.5498E  + 

00 

O.OOOOE +00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

-1.5498E  +  00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

-l.OOOOE-04 

O.OOOOE +00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE +  00 

O.OOOOE +00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-l.OOOOE-  04 

O.OOOOE  + 

00 

O.OOOOE  + 

00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

O.OOOOE  +  00 

-l.OOOOE- 

04 

Bi  = 


0.0000E+  00 
0.0000E+  00 
0.0000E+  00 
0.0000E+  00 
O.OOOOE+00 
0.0000E+  00 
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